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Abstract 

In this paper, we consider systems of algebraic and non-linear partial differential equations and 
inequations. We decompose these systems into so-called simple subsystems and thereby partition 
the set of solutions. For algebraic systems, simplicity means triangularity, square- freeness and non- 
vanishing initials. Differential simplicity extends algebraic simplicity with involutivity. We build 
upon the constructive ideas of J. M. Thomas and develop them into a new algorithm for disjoint 
decomposition. The given paper is a revised version of Bachler et al. (2010) and includes the 
proofs of correctness and termination of our decomposition algorithm. In addition, we illustrate 
the algorithm with further instructive examples and describe its Maple implementation together 
with an experimental comparison to some other triangular decomposition algorithms. 

Keywords: disjoint triangular decomposition, simple systems, polynomial systems, differential 
systems, involutivity 



1. Introduction 

Nowadays, triangular decomposition algorithms, which go back to the characteristic set method 
by Ritt (1950) and Wu (2000), have become powerful tools for investigating and solving systems 
of multivariate polynomial equations. In many cases these methods are computationally more 
efficient than those based on construction of Grobner bases. For an overview over triangular 
decomposition methods for polynomial and differential-polynomial systems we refer to the tutorial 
papers by Hubert (2003a,b) and to the bibliographical references therein. 

Among numerous triangular decompositions the Thomas one stands by itself. It was sug- 
gested by the American mathematician Thomas (1937, 1962) and decomposes a finite system 
of polynomial equations and/or inequations into finitely many triangular subsystems, which he 
called simple. The THOMAS decomposition splits a given quasi-affine variety into a finite number 
of quasi-affine varieties defined by simple systems. Unlike other decomposition algorithms, the 
Thomas decomposition always yields a disjoint decomposition of the solution set. 

Wang was the first to design and implement an algorithm that constructs the THOMAS decom- 
position (cf. Wang (1998, 2001); Li and Wang (1999)). For polynomial systems he implemented his 
algorithm in Maple (cf. Wang (2004)) as part of the software package epsilon (cf. Wang (2003)), 
which also contains implementations of a number of other triangular decomposition algorithms. 
Delliere (2000) has shown that the "dynamic constructible closure" introduced in the thesis by 
Gomez Diaz (1994) can be modeled using simple systems. Nonetheless, according to the remark 
after (Delliere, 2000, Thm. 5.2), simple systems are more general. 

Every simple system is a regular system and its equations form a regular chain. The Regular- 
Chains package (cf. Lemaire et al. (2005)) includes procedures to decompose the solution set of the 
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input by means of regular chains (if the input only consists of equations) or regular systems. How- 
ever, the Thomas decomposition differs noticeably from this decomposition, since the Thomas 
decomposition is finer and demands disjointness of the solution set. For a detailed description of 
algorithms related to regular chains, we refer the reader to Moreno Maza (1999). 

The disjointness of the Thomas decomposition combined with the structural properties of 
simple systems provide a useful platform for counting solutions of polynomial systems. In fact, the 
THOMAS decomposition is the only known method to compute the counting polynomial introduced 
by Plesken (2009a). We refer to §2.3 for details on this structure, counting and their applications. 

During his research on triangular decomposition, Thomas was motivated by the Riquier- 
Janet theory (cf. Riquier (1910); Janet (1929)), extending it to non-linear systems of partial 
differential equations. For this purpose he developed a theory of (THOMAS) monomials, which 
generate an involutive monomial division nowadays called Thomas division (cf. Gerdt and Blinkov 
(1998a)). He gave a recipe for decomposing a non-linear differential system into algebraically 
simple and passive subsystems (cf. Thomas (1937)). A modified version of the differential Thomas 
decomposition was considered by Gerdt (2008) with its link to the theory of involutive bases (cf. 
Gerdt and Blinkov (1998a); Gerdt (2005, 1999); Seiler (2010)). In this decomposition, the output 
systems are jANET-involutive in accordance to the involutivity criterion from Gerdt (2008) and 
hence they are coherent. For a linear differential system it is a Janet basis of the corresponding 
differential ideal, as computed by the Maple package Janet (cf. Blinkov et al. (2003)). 

The differential THOMAS decomposition differs from that computed by the Rosenfeld-Grob- 
NER algorithm (cf. Boulier et al. (2009, 1995)). The latter decomposition forms a basis of the dif- 
falg, DifferentialAlgebra and BLAD packages (cf. Boulier and Hubert (1996-2004); Boulier (2004-2009)). 
Experimentally, we found that these three packages are optimized and well-suited for ordinary 
differential equations. Furthermore, epsilon also allows to treat ordinary differential systems. 
Bouziane et al. (2001) mentions another implementation not available to the authors. However, 
all these methods give a zero decomposition, which, unlike the Thomas decomposition, is not 
necessarily disjoint. 

In the given paper we present a new algorithmic version of the Thomas decomposition for 
polynomial and (ordinary and partial) differential systems. In this unified algorithm, only two 
changes to the algebraic version are necessary to adapt it for the treatment of differential systems. 
We briefly describe our implementation of this algorithm in Maple. 

This paper is organized as follows. In §2, we present the algebraic part of our algorithm for 
the Thomas decomposition with its main objects defined in §2.1. In §2.2, we describe the main 
algorithm and its subalgorithms, followed by the correctness and termination proof. Decomposition 
of differential systems is considered in §3. Here, we briefly introduce some basic notions and 
concepts from differential algebra (§3.1) and from the theory of involutive bases specific to Janet 
division (§3.2). In §3.3, we present our version of the differential pseudo reduction and, building 
upon it, the definition of differential simple systems. Subsection §3.4 contains a description of the 
differential Thomas decomposition algorithm and the proof of its correctness and termination. 
Some implementation issues are discussed in §4, followed by a comparison of our implementation 
to some other implementations of triangular decompositions with the help of benchmarks. 

2. Algebraic Systems 

This section introduces the concepts of simple systems and the Thomas decomposition for 
algebraic systems. These concepts are based on properties of the set of solutions of a system. We 
conclude the section with an algorithm for constructing a Thomas decomposition. 

Example 2.1. We give an easy example of a Thomas decomposition. Consider the equation 

P = x 3 + (3y + l)x 2 + (3y 2 + 2y)x + y 3 = . 
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A THOMAS decomposition of {p = 0} is given by: 

51 :={ x 3 + (3y + l)x 2 + (3y 2 + 2y)x + y 3 = 0, 

27y 3 - 4y ? 0} 

5 2 := { 6x 2 + {-27y 2 + I2y + 6)x - 3y 2 + 2y = 0, 

27y 3 - Ay = 0} 



The picture shows the solutions of {p = 0} in the real affine plane. The cardinality of the 
fibers of the projection onto the y-component depends on y. However, if we consider all solutions 
in the complex plane, this cardinality is constant within each system, i.e., 3 and 2 in S\ and S2, 
respectively. This property is formalized in the definition of simple systems. 




2.1. Preliminaries 

Let Fbea computable field of characteristic and R := F[x\, . . . , x n ] be the polynomial ring 
in n variables. A total order < on {l,xi, . . . ,x n } with 1 < Xj for all i is called a ranking. The 
indeterminate x is called leader of p G R if x is the <-largest variable occurring in p. 2 In this 
case we write ld(p) = x. If p G F, we define ld(p) = 1. The degree of p in ld(p) is called main 
degree of p (mdeg(p)) and the leading coefficient init(p) £ F[ y \ y < ld(p) ] of ld(p) mdeg ^ p ^ in p 
is called initial of p. 

For a G F , where F denotes the algebraic closure of F, define the following evaluation 
homomorphisms: 

a : F[xi , . . . , x n ] -> F : xi i-)- a-i 



" " i,:/ '- r; ■''■ n : { ZZZ otherwise 

Given a polynomial p G R, the symbols p = and p^ denote the equation p = and inequation 
p / 0, respectively. A finite set of equations and inequations is called an (algebraic) system over 
R. Abusing notation, we sometimes treat p = or p^ as the underlying polynomial p. A solution 
of p = or p^ is a tuple a G f" with (f> a (p) — or (f> a (p) / 0, respectively. We call a G f" a solution 
of a system S, if it is a solution of each element in S. The set of all solutions of S is denoted by 
&o^S). 

The subsets of all equations p = G S and all inequations G S are denoted by S = and S^, 
respectively. Define S x :— {p € S \ ld(p) = x}. In a situation where it is clear that | S x I — 1 7 we 
also write S x to denote the unique element of S x . The subset S <x := {p G S* | ld(p) < x} is a 
system over F[y \ y < x ]. 

The Thomas approach uses the homomorphisms 4>< x ^ to treat each polynomial p G S x as the 
family of univariate polynomials 0<a:,a(p) G -^t 21 ] f° r a G 6ol(/S , < x ). This idea forms the basis of 
our central object, the simple system: 

Definition 2.2 (Simple Systems). Let S be a system. 

1. S is triangular if \S Xi < 1 V 1 < i < n and S n {c = , | c G F} = 0. 

2. S has non-vanishing initials i/0 a (init(p)) ^ V a G (3o[(S f < a : i ) anc?p G S Xi for 1 < i < n. 

3. S is square-free 3 if the univariate polynomial 4>< Xil a(p) G F[xi] is square-free\/& G &ol(S <Xi ) 
and p G S Xi for 1 < i < n. 



2 In the context of triangular decompositions, the leader is usually called main variable. The term leader is 
used in Thomas (1937) and has later been adopted in differential algebra. 

3 Square-freeness has an important side-effect in the differential case. A square-free polynomial and its separant 
have no common roots. Thus, the separants do not vanish on solutions of the lower-ranking subsystems. 
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4-- S is called simple if it is triangular, has non-vanishing initials and is square-free. 



Properties (2) and (3) are characterized via solutions of lower-ranking equations and inequa- 
tions. However, the Thomas decomposition algorithm does not calculate roots of polynomials. 
Instead, it uses polynomial equations and inequations to partition the set of solutions of the 
lower-ranking system to ensure the above properties. 

Remark 2.3. Every simple system has a solution. In particular, ifh€ &ol(S <x ) and S x is not 
empty, then 4><x,h(S x ) is a univariate polynomial with exactly rndeg(S x ) distinct roots. When S x 
is an equation, each solution b £ &ol(S <x ) extends to a solution (b, a) 6 &ol(S< x ) with mdeg(S x ) 
possible choices a € F. Otherwise, all but finitely many a € F yield a solution (b,a) G &ol(S< x ), 
because an inequation S x excludes mdeg(S , x ) different a and S x = imposes no restriction on a. 

Conversely, if (a±, . . . , a n ) G ©ol(S) where S is a system over F[x\, . . . , x n ] with X\ < . . . < x n , 
then (ax, . . . , aj) e &ol(S< Xi ). 

To transform a system into a simple system, it is in general necessary to partition the set of 
solutions. This leads to a so-called decomposition into simple systems. 

Definition 2.4. A family (Si) 7 ^ is called decomposition of S if &ol(S) — U"=i SoI(Si). A 
decomposition is called disjoint if 6ol(Si) H &ol(Sj) = V i ^ j. A disjoint decomposition of a 
system into simple systems is called (algebraic) THOMAS decomposition. 

For any algebraic system S, there exists a Thomas decomposition (cf. Thomas (1937, 1962), 
Wang (1998)). The algorithm presented in the following section provides another proof of this 
fact. 

Example 2.5. We compute a Thomas decomposition of { (p :— ax 2 + bx + c) _} C Q[a, b, c, x] 
with respect toa<b<c<x. We highlight the highest power of the leader by underlining it. 

First, we ensure that the initial init(p) of p is not zero. Therefore, we insert (init^))^ = (a)_^ 
into the system. Since we restricted the solution set of this system, we also have to consider 
the system {p=,(<l) = }, which simplifies to {(bx + c)_ , (a) = }. Similarly, we add (b)^ to ensure 
init(6a; + c) ^ and get the special case system {(c)_ , (b)_ , (a) = }. Up to this point, we have three 
systems, where the second and third one are easily checked to be simple: 

x | (ax 2 + bx + c) = x. (bx + c) = • 



c ? 



b " (b) 



(s)= 

(b)-_ 



a I (a)_£ a I (a) = ! 

Second, we ensure that p is square-free by insertion of (4ac — £> 2 ) , into the first system. Again, 

we also need to consider the system {(p) = , (4ac — b 2 )_ , (a)^}- As p is a square in this system, 
we can replace it by its square-free part 2ax — b. Now, all systems are easily verified to be simple 
and we obtain the following THOMAS decomposition: 



x . (ax z 



bx 



Uac-b 2 



x | (2ax - b) = 
(Aac~b 2 ) = 



X m (bx + C)__ 



c <> 



a I (a) 



b o ® # 

a i (a) = 



(s)= 

(b)-_ 
a I (a)__ 
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2.2. Algebraic Thomas Decomposition 

This section presents our main algorithm for algebraic systems and its subalgorithms. The 
algorithm represents each system as a pair consisting of a candidate simple system and a queue 
of unprocessed equations and inequations. During each step, the algorithm chooses a suitable 
polynomial from the queue, pseudo-reduces it and afterwards combines it with the polynomial 
from the candidate simple system having the same leader. In this process, the algorithm may 
split the system, i.e., add a new polynomial into the queue as an inequation and at the same time 
create a new subsystem with the same polynomial added to the queue as an equation. This way, 
we ensure that no solutions are lost and the solution sets are disjoint. The algorithm considers 
a system inconsistent and discards it when an equation of the form c— with c G F \ {0} or the 
inequation 0^ is produced. 

We consider a system S as a pair of sets (St, Sq), where St represents the candidate simple 
system and Sq is the queue. We require St to be triangular and thus (St)x denotes the unique 
equation or inequation of leader x in St, if any. Moreover, St must fulfill a weaker form of the 
other two simplicity conditions, in particular, in conditions (2.2)(2) and (3), the tuple a can be a 
solution of {St)<x U {Sq)< x instead of just {St)<x- Obviously, Sq = implies simplicity of S. 

From now on, let prem be a pseudo remainder algorithm 5 in R and pquo the corresponding 
pseudo quotient algorithm. To be precise, if p, q G R with ld(p) = \d(q) = x, then 

m ■ p = pquo(p, q, x) ■ q + prem(p, q, x) (1) 

holds, where deg x {q) > deg x (prem(p, q, x)), ld(m) < x and m | init(q) fc for some k G Z>o- Note 
that </> a (init(p)) ^ and </> a (init(g)) ^ imply a (pquo(p, q, x)) ^ and </> a (m) ^ 0. 
The following algorithm employs prem to reduce a polynomial modulo St- 

Algorithm 2.6 (Reduce). 

Input: A system S, a polynomial p G R 

Output: A polynomial q with <^ a (p) = if and only if <fi a (q) — for each a G &ol(S). 
Algorithm: 
1: x <— ld(p); q <— p 

2: while x > 1 and {St) x is an equation and mdeg((j) > mdcg((Sr)a;) do 
3: q <- prem(<3[, (S T ) X ,%) 
4: x <r- ld(q) 
5: end while 

6: if x > 1 and Reduce^, init(q)) = then 
7: return Reduce^, q - init(g)^ mdeg ^)) 
8: else 

9: return q 
10: end if 

Proof (Correctness). There exist m G R \ {0} with ld(m) < ld(p) and <p a (m) ^ for all 
a G ©o[(5<id( p )) such that 

Reduce(5,p) = mp — c y ■ {St) v 

y<ld(p) 

with c y G R and ld(cj,) < ld(p) if {St) v is an equation and c y — otherwise. This implies 
a (Reduce(S») = <p a (m) a (p) - ^2 ^ a ( c f) M{St) v ) 



4 This approach has been adapted from Gerdt and Blinkov (1998b), where T was an intermediate Janet basis 
and Q a queue of new prolongations to be checked. A similar approach was later used for triangular decompositions 
in Moreno Maza (1999). 

5 In our context prem does not necessarily have to be the classical pseudo remainder, but any sparse pseudo 
remainder with property (1) will suffice. 
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and therefore 4> a (p) = if and only if </> a ( Reduce (5, p)) = 0. 



□ 



Note that this algorithm only uses the equation part of the triangular system in S, i.e. ■ 
For ease of notation in the following algorithms, we write Reduce(5,p) instead of Reduce(S^ , p) . 

A polynomial p reduces to q modulo St if Reduce(5, p) = q. A polynomial is reduced 
modulo St if it reduces to itself. 

The Reduce algorithm differs slightly from the classical prem(p, S^f) as defined in Aubry et al. 
(1999). While prem(p, S^ ) fully reduces p modulo all variables, Reduce(S',p) only reduces modulo 
the leader and ensures that the initial of the reduced form doesn't vanish. Performing Reduce(S*, p) 
in combination with a full coefficient reduction (see also §4.2) is the same as computing prem(p, S^)- 
It is therefore possible to replace Reduce(5, p) with prem(p, S^) in the following algorithms. Our 
approach adds some flexibility, as we can choose to omit a full reduction in an implementation. 
In particular, if a polynomial does not reduce to zero, we can determine that without performing 
a full prem reduction. We apply this multiple times in our implementation, most prominently in 
Algorithm (2.18). However, if a polynomial reduces to zero, Reduce has no advantage over prem. 

Later, we will use the following facts about the Reduce algorithm. 

Remark 2.7. Let q = Reduce(5,p) ^ 0. 

1. //S'id(g) is an equation, then vadeg{q) < mdeg(S\d( q )) ■ 

2. Reduce(S,init(Reduce(S>))) ^ 0. 

3. ld(g) < \d{p) and if\d(q) = \d(p), then mdeg(<7) < mdeg(p). 

The result of the Reduce algorithm does not need to be a canonical normal form, however, the 
algorithm recognizes polynomials that vanish on all solutions: 

Corollary 2.8. Let p G R with ld(p) = x. Reduce(5,p) = implies <fi a (p) = V a E ©ol(5'< x ). 

Proof. For all a <E &ol(S< x ), it holds that <f> a {p) = if and only if (f> a (Reduce(S,p)) = 0. The 
statement follows from </> a ( Reduce^, p)) = 4> a (0) = 0. □ 

The converse of this corollary doesn't hold in general. Thus, we provide two weaker statements 
in the following remark. 

Remark 2.9. Let p and x as in Corollary (2.8). 

1. If {Sq)< x = 0, i.e., S< x — (St)<x is simple, then Reduce(5,p) =/= implies 3 a G &ol(S< x ) 
such that 4> a (p) =/= 0. 

2. If (S Q )= X = and Reduce(S,p) ^ hold, then either &ol(S <x ) = or 3 a e 6ol{S <x U 
{{St) x }) such that 4> a (p) ^ 0. 

Proof. We only prove the second part, as the first part easily follows. 

Let (S Q )= X = 0, Reduce(S» ^ and |Sol(5< K )| > 0. First, as \d((S T ) x ) ^x and 
mdeg((ST) x ) > 0, for each a 6 &ol(S <x ), the univariate polynomial <fi< x , a ((ST) x ) G F[x] has 
positive degree. Thus \&ol(S <x U {{S T )x})\ > 0. 

Let <fi a (p) = V a G &ol(S <x U {{St) x }) (*)■ Then (St) x is an equation and deg^p) > 
deg x ((ST) x ) and therefore p ^ Reduce(5,p). In fact, (*) further implies ld(Reduce(S', p)) < x, as 
otherwise deg :E (Reduce(S', p)) > deg x ((ST) x ) would hold. By repeating the previous arguments, 
we can inductively conclude ld(Reduce(S', p)) = 1. As 4> a (p) = 0, we conclude Reduce(S*, p) = 0, a 
contradiction. □ 

The first part of this remark in conjunction with Corollary (2.8) implies (Wang, 1998, Thm. 4). 
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Example 2.10. Reduce qi := x 2 + y 2 x + x + y modulo the simple system on the left. 
S x = (yx 2 -1) = x 2 + y 2 x + x + y =: qi 

J yqi-S x 

{y 3 + y)x + y 2 + 1 =:q 2 



yiSy = (y 2 + l) = ^y 2 + l=:q 3i _> y J + y = init(q 2 ) 

/mit(g 2 ) - y-S y 

0^ 

In the first reduction step, q\ is pseudo-reduced modulo S x . The result q 2 still has leader x, but 
a main degree smaller than S x . We determine that the initial of qi reduces to and remove the 
highest power of x from q 2 . The resulting polynomial q 3 now pseudo-reduces to modulo S y , i.e. 
Reduce({S x ,S y }, qi ) = 0. 

Now, we examine all splitting methods needed during the algorithm. We will use the following 
one-liner as subalgorithm for the splitting subalgorithms. 

Algorithm 2.11 (Split). Input: A system S, a polynomial p G R 
Output: The disjoint decomposition {S U {p^} , S U {p=}) of S . 
Algorithm: 

i: return {{St, Sq U {p*}) , {St, Sq U {p=})) 

For a better understanding of the following splitting subalgorithms we first need to explain 
how they are applied in the main algorithm. Each step of the algorithm treats a system S as 
follows. An equation or inequation q is chosen and removed from the queue Sq. Then we reduce 
q modulo St- For the simplicity properties to hold w.r.t. q it is necessary to add inequations 
to S. To accomplish this, we pass S together with q to the splitting subalgorithms. Each such 
subalgorithm returns two systems. The first system Si contains an additional inequation. The 
second system S 2 contains a complementary equation, q is added back into the queue of S 2 , and 
S 2 is put aside for later treatment. In each case {Si U {q}, S 2 ) is a disjoint decomposition of the 
original system S U {q}. Then Si and q may be subjected to further splitting algorithms and 
eventually q is added into the candidate simple system. 

The first splitting algorithm we consider is InitSplit, which is concerned with property (2.2)(2). 

Algorithm 2.12 (InitSplit). Input: A system S, an equation or inequation q with ld{q) = x. 
Output: Two systems Si and S 2 , where {Si U {q}, S 2 ) is a disjoint decomposition of S U {q}. 
Moreover, <^ a (init(#)) 7^ holds for all a G ©o[(Si) and a (init(g)) = for all a G &ol{S 2 ). 
Algorithm: 

i: (5i,5 2 ) <-Split(S,init( 9 )) 
{S 2 ) Q <- {S 2 ) Q U {q} 
return {Si,S 2 ) 

For the further splitting algorithms, we need some preparation. In Definition (2.2) we consider 
a multivariate polynomial p as the family of univariate polynomials 0<id(p),a(?O- For ensuring 
triangularity and square-freeness, we have to compute the gcd (greatest common divisor) of two 
polynomials, which in general depends on a. Subresultants provide a generalization of the Eu- 
CLiDean algorithm and enable us to take the tuple a into account. 

Definition 2.13. Let p,q G R with ld(p) = ld(g) = x, deg x {p) — d p > deg x {q) — d q . We 
denote by PRS(p, q, x) the subresultant polynomial remainder sequence (see Habicht (1948), 
(Mishra, 1993, Chap. 7), (Yap, 2000, Chap. 3)) of p and q w.r.t. x, and by PRSi(p, q, x), i < 
d q the regular polynomial of degree i in PRS(p, q, x) if it exists, or otherwise. Furthermore, 
PRS dp (p, q, x) := p, PRS dg {p, q, x) := q and PRSi{p,q,x) := 0, d q < i < d p . 
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Define teSi(p, q, x) := init (PRS; (p, q, x)) for < i < d p , tesd p {p, q, x) := 1 and reso(p, q, x) := 
PRSo (p, (?, a;) . 7Vo£e thatreso(p,q,x) is the usual resultant. 6 

The initials of the subresultants provide conditions to determine the degrees of all possible gcds. 
Using these conditions, we describe the splittings necessary to determine degrees of polynomials 
within one system. 

Definition 2.14. Let S be a system and p%,P2 € R with \d(pi) = ld(p2) = x. If |©oI(5< x )| > 0, 
we call 

i := min{i G Z> | 3 a G &ol(S <x ) such that deg x (gcd{(f> <x . a (pi),(j) <x . a {p2))) = i} 

the fiber cardinality of pi and p2 w.r.t. S. Moreover, if {Sq)^ x — 0, then 

i' := min{i G Z> | Reduce(reSj(pi,p2, x), St) — V j < i and Reduce(resj(f>i,p2i x), St) ^ 0} 

is the quasi fiber cardinality of p\ and p2 w.r.t. S. A disjoint decomposition (Si,^) of S such 
that 

1- deg a ,(gcd(0 <2; , a (pi),(/> <2; , a (p2))) = i V a G &ol ((Si) <x ) 

2. deg 2; (gcd(0 <2; , a (pi),</> <2; , a (p2))) >!Vae6ol((S 2 y 

is called i-th fibration split of pi and p2 w.r.t. S. A polynomial r G R with ld(r) = x such that 
deg x (r) — i and 

(f><x.a( r ) ~ gcd(0 <a; , a (pi), (/> <x ,a.(p2)) V a G 6o^ ((Si) <x ) 
is called i-th conditional greatest common divisor of pi and p 2 w.r.t. S, where p ~ q if and 
only if p G (F \ {0})q. Furthermore, q G R with ld(q) = x and deg x (q) = deg x (pi) — i such that 

~ au ( ^^ V a e eo[ ((^i)<,) 

gcd(0 <a; , a (pi), <p <x , a (P2)) 

is called the i-th conditional quotient of pi by P2 w.r.t. S. By replacing 4 , <x,a.(P2) i n the above 
definition with ^(0<x,a(pi)); we get an i-th square-free split and i-th conditional square-free 
part of pi w.r.t. S. 

Example 2.15. Consider the system S :— {(x 3 + y) = } and the polynomial q :— x 2 + x + y+l with 
y < x. Compute reso(S x , q) = y 3 + 7y 2 + 5y + 1, res\{S x ,q) = —y and res2(S x ,q) = 1. The fiber 
cardinality is of S x and q w.r.t. S isO. A zeroth fibration split is given by Si :— SU{(reso(S x , 
and S2 '■= S U {(reso(S x , q))=}. The fiber cardinality w.r.t. S2 is 1. A first fibration split is given 
by 52,1 := S2U{(— y)^} and S 2, 2 '■= 5 U {(—?/) = }. Note in this case that &o\(S 2,1) = &oi{Si) and 
6o[(S f 2,2) = %■ A zeroth conditional quotient of S x and q is S x . A first conditional gcd and first 
conditional quotient are —yx + 2y + 1 and y 2 x 2 + (fly 2 + y)x + Ay 2 + Ay + 1, respectively. 

It is in general hardly possible to compute the fiber cardinality directly. However, in the case 
where the quasi fiber cardinality is strictly smaller than the fiber cardinality, the corresponding 
fibration split will lead to one inconsistent system, and one where the quasi fiber cardinality is 
increased. 

Lemma 2.16. Let \<Sol(S< x )\ > and (Sq)^ x — 0. For pi, P2 as in Definition (2.14) with 
a (init(pi)) 7^ V a G &ol(S <x ) and mdeg(pi) > mdeg(p2) 7 let i be the fiber cardinality of pi and 
P2 w.r.t. S and i' the corresponding quasi fiber cardinality. Then 

1 < 1 

where the equality holds if and only if \&ol(S^ x U {resi>(pi,p2,x)^})\ > 0. 

6 These definitions are slightly different from the ones cited in the literature ((Mishra, 1993, Chap. 7), (Yap, 
2000, Chap. 3)), since we only use the regular subresultants. However, it is easy to see that all theorems from 
(Mishra, 1993, Chap. 7) we refer to still hold for i < d q . 
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PROOF. Let a e &ol(S <x ), mdeg(pi) > mdeg(p 2 ), d Pl := deg x (p 1 ) = deg x (0 <x , a (pi)), d P2 := 
deg x (p2) and d P2t8L := deg x (</><;r !a (p2))- If i < max(d Pl , <i P2 , a ) — 1 = d Pl — 1, then (Mishra, 1993, 
Thm. 7.8.1) implies 

<A<x,a(PRSi(pi,P2,2c)) ~ PRS l (0 <a;ia (pi),0 <:C:a (p2),a;) (2) 

and 

</> a (reSi(pi,p 2 ,a;)) = reSi(0 <X:a (pi),0< x , a (p 2 ),a;) = . (3) 

Conditions (2) and (3) by definition also hold for the trivial cases d P2 < i < d Pl . 

For all indices j < i', Corollary (2.8) and the fact Reduce(reSj ;(pi,P2> x), St) = imply 
(f) a (resj(pi,p2,x)) = 0. By (2) and (3), reSj (0<x, a (Pi))<^<a;,a(P2),2c) = follows. We apply 
(Mishra, 1993, Thm. 7.10.5) successively and get PRJ3,(^>< Xia (pi), (j>< x ,&{P2), x) = 0. Thus, 

deg x (gcd(0 <x , a (>i),0 <x , a (p2))) > i' (4) 

holds. This implies i' < i. 

Equality in (4) holds if and only if there exists a 6 &ol(S <x ) such that a (reSi' (pi,P2, x)) ^ 0. 
Therefore, i = %' if and only if \&ol(S <x U {resi/ (pi,P2, a0^})| > 0. □ 

The above lemma doesn't apply if both polynomials have the same degree. In this case, both 
polynomials must have non-vanishing initials, as shown in the following corollary. 

Corollary 2.17. Let \&ol(S <x )\ > and (£q)< x — 0. For polynomials p\, p2 as in Definition 
(2.14) with a (init(pi)) ^ and a (init(j>2)) 7^ V a 6 &ol(S <x ), let i be the fiber cardinality of 
pi and pi w.r.t. S and i' the quasi fiber cardinality of pi and prem(p2,Pi, x) w.r.t. S. Then 

i' < i 

with equality if and only if \&ol (S <x U {resj'(pi, prem(p2>Pi) x), x)^;})\ > 0. 

Proof. Let a 6 &ol(S <x ). By the assumption on the initials, (Mishra, 1993, Corr. 7.5.6) implies 
0<x,a(prem(p 2 ,Pi,aO) = prem(</>< x , a (p 2 ), <£<s,a(Pi), x). The univariate polynomials 4><x, a (j>i) and 
4><x,&(jp2) have the same gcd as <p< Xl a(pi) and prem(</> <:Cja (p2), <fr<x,a(pi), x). We can therefore 
replace P2 with prem(p2,Pi, x ) in Lemma (2.16). □ 

The following algorithm computes the quasi fiber cardinality of two polynomials. 

Algorithm 2.18 (ResSplit). Input: A system S with {Sq)^ x — 0, two polynomials p,q £ R with 

ld(p) = ld(q) = x, mdeg(p) > mdeg(q) and </> a (init(p)) ^ for all a e &ol(S <x ). 

Output: The quasi fiber cardinality i of p and q w.r.t. S and an i-th fibration split (SijS^) of p 

and q w.r.t. S . 

Algorithm: 

1: i «— min{i € Z> | Reduce(reSj(p, q, x), St) = V j < i and Reduce(reSi(p, q, x),St) ^ 0} 
2: return (i, Si, S2) := (i, Split(5,reSj(p, g, a;))) 

Proof (Correctness). Assume \&ol((Si) <x )\ > 0, 1 = 1,2, as the statement is trivial otherwise. 

Let a S &ol((Si) <x ). The polynomial g := PKSi(4> <x ^ SL (p) , (f> <x ^ a (q) , x) is not identically zero, 
due to (init(p))^ = (resj(p, <Z, x))^ € {Si)q. The degree of g is i and g ~ gcd(0 <a;ia (p), (f> <x<gL (q)), 
as discussed in the proof of Lemma (2.16). 

Let a € 6o[((S , 2 )< K ). (Mishra, 1993, Thm. 7.10.5) and (init( 5 )) = = (reSi(p, q, x))= e {S 2 )q 
imply 3 = 0. Therefore, deg x (gcd(<£ <x , a (p), 4><x,a(<l))) > i- □ 

We apply the fiber cardinality and fibration split to compute a greatest common divisor of an 
existing polynomial in St and another polynomial. 
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Algorithm 2.19 (ResSplitGCD). Input: A system S with (Sq)^ x — 0, where (St)x * s an equa- 
tion, and an equation q— with \d(q) = x. Furthermore, mdeg(<7) < mdeg((S , 7') x ). 
Output: Two systems S\ and S 2 and an equation q— such that: 

a) S 2 = S 2 U {q} where (s±, S2J is an i-th fibration split of (St)x and q w.r.t. S , 

b) q is an i-th conditional gcd of (St)x and q w.r.t. S, 

where i is the quasi fiber cardinality of p and q w.r.t. S. 
Algorithm: 

i: (i,S u S 2 ) <- ResSplit(S, (S T ) x ,q) 

2: (S 2 ) Q <- (S 2 ) Q U {q} 

3: return S%, S 2 , PRSi((S T )x, q, x)= 

Proof (Correctness). Property a) follows from Algorithm (2.18) and line 2. Property b) was 
already shown in the correctness proof of Algorithm (2.18). □ 

Note that i > is required in this i = would yield an inconsistency. Therefore, before 

calling ResSplitGCD. we will always ensure this condition in the main algorithm by incorporating 
the resultant of two equations into the system. 

The following algorithm is similar. But instead of the gcd, it returns the first input polynomial 
divided by the gcd. 

Algorithm 2.20 (ResSplitDivide). Input: A system S with (Sq)< x = and two polynomials p, 
q with ld(p) = ld(q) — x and a (init(p)) ^ for all a E &ol(S <x ). Furthermore, if mdeg(p) < 
mdeg(g) then (f> a (mit(q)) ^ 0. 

Output: Two systems Si and S 2 and a polynomial p such that: 

a) S 2 = S 2 U {q} where is an i-th fibration split of p and q' w.r.t. S, 

b) p is an i-th conditional quotient of p by q' w.r.t. S, 

where i is the quasi fiber cardinality of p and q' w.r.t. S , with q' = q for mdeg(p) > mdeg(g) and 

q' = prem(q,p, x) otherwise. 

Algorithm: 

1: i/mdeg(p) < mdeg(g) then 

2: return ResSplitDivide(5 f , p, prem(q,p, x)) 

3: else 

4.- (i,Si,S 2 ) <- ResSpWt (S,p,q) 
5: if i > then 

6: p ^— pquo(p, PRSi(p, prem(g,p, x), x), x) 
7: else 

8: p <— p 

9: end if 

10: [S 2 )q <- (S 2 ) Q U {q} 

11: return S\,S 2 ,p 
12: end if 

Proof (Correctness). According to Corollary (2.17), we can without loss of generality assume 
mdeg(p) > mdeg(g). 

Property a) follows from Algorithm (2.18) and line 10. For all a £ ©oI(5i), the following holds: 
If i = 0, then deg x (gcd((f) <Xya (p),(j) <x , a (q'))) = and thus 4>< x ,a.{p) shares no roots with (j)< x ^(q'). 
Now let i > 0. Formula (1) implies 

m-p — p- PRSi (p, q', x) + prem (p, PRSi [p, q', x) , x) . 
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Due to (Mishra, 1993, Cor. 7.5.6) and (2), (3) there exist ki,k 2 G F\ {0} such that 
= <P<xAp) ■ ^<a;,a (PRS 4 (p, a;)) + <a ., a (prem(p, PRSi (p, <j, x) , a;)) 

= 0<^,a(p) • fci PRSi (</> <a: ,a(p), <t><x,a.{q), x) + k 2 prem((j) <x , a (p) , PRS; (0< x , a (p)i 0<o:,a(?)> x), x) 



divides 0< x ,a(p) 

= 0<x,a(p) ' ^1 gcd(^ <x>a (p), 0< a ,a(<?)) + . 

Thus, we obtain property b) from 



*(P) 



gcd(0< X:a (p),(/)< 2: , a (g)) 

and deg x (0 <x , a (p)) = deg a ,(0 <x ^ a (p)) - deg x (gcd(^ <x>a (p), 0<x,a(?))) = deg x (p) - i. □ 

Applying the last algorithm to a polynomial p and -q^j-^P yields an algorithm to make p 
square-free. We present it separately for better readability of the main algorithm. 

Algorithm 2.21 (ResSplitSquareFree). Input: A system S with (Sq)< x = and a polynomial p 
with ld(p) = x and a (init(p)) ^ for all a £ &ol(S <x ). 
Output: Two systems Si and S 2 and a polynomial r such that: 

a) S2 = S2 U {p} where (Si,S%\ is an i-th square-free split of p w.r.t. S, 

b) r is an i-th conditional square-free part of p w.r.t. S, 

where i is the quasi fiber cardinality of p and -^p w.r.t. S . 
Algorithm: 

v. (i,S u S 2 ) <- ResSplit(S,p,,§p) 

2: if i > then 

3: r <- pquo (p, PRSi (p, J^p, x) , x) 

4: else 
5: r <— p 

6: end if 

7: (S 2 ) Q <- (S 2 ) Q U {p} 

S: return Si,S2,r 



Proof (Correctness). Since 4> <x a (-§-p) — -§-<t><x a(p), an i-th square-free split of p is an i-th 
d_ 

y.r-< 



fibration split of p and -J^-p. The rest follows from the proof of Algorithm (2.20). □ 



In all ResSplit-based algorithms, (Sq)< x = is required. This ensures that all equations of a 
smaller leader than x will be respected by reduction modulo St- The order in which polynomials 
are treated by the main algorithm must therefore be restricted. 

Definition 2.22 (Select). Let Pfi nite (M) be the set of all finite subsets of a set M. A selection 
strategy is a map 

Select : Pfi n ite({p=,P^ I P G R}) — ► {p=,P^ I P G R} : 

Q 1 — > qeQ 

with the following properties: 

1. //Select(<5) = q— is an equation, then Q^j^ = 0- 

2. //Select(<5) = is an inequation, then Q<w ? ) = fi- 
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We demonstrate that these conditions are necessary for termination of our approach, by giving 
an example where we violate them. 

Example 2.23. Consider R := F[a, x] with a < x and the system S with St ■— and Sq := 
{(x 2 — «)=}■ To insert (x 2 — a) = into St, we need to apply the ResSplitSquareFree algorithm: We 
calculate reSo(a; 2 — a, 2x, x) = —4a, reSi(:c — a, 2x, x) = 2 and ies 2 (x 2 — a, 2x, x) = 1 according to 
Definition (2.13). The quasi fiber cardinality is and we get the two new systems S\, S 2 with 

(5i) T = {(x 2 - o)=},(Si) = {(-4a) # } and (S 2 ) T = 0, (S 2 )q = {(x 2 - o)=, (-4o)=} . 

We now consider what happens with S 2 : If we select (x 2 — a)— as the next equation to be 
treated, in violation of the properties in Definition (2.22), ResSplitSquareFree will split up S 2 into 
S21, S22 with 

(S 2 ,i)t = {(x 2 - a)=}, (S 2A ) Q = {(-4a) #! (-4o)=} 

and 

(S 2 , 2 ) T = 0, (S 2 ,2)q = {(x 2 - a)=, (-4a) = , (-4o)=} . 
j4s 1S2 = S 2j 2, this will lead to an endless loop. 

The following trivial algorithm inserts a new equation into St- It will be replaced with a 
different algorithm in §3 when the differential THOMAS decomposition is considered. 

Algorithm 2.24 (InsertEquation). Input: A system S and an equation r— with ld(r) = x satisfy- 
ing </> a (init(r)) ^ and <p< x .a.{r) square-free for all a G &ol(S <x ). 
Output: A system S where r = is inserted into St- 
Algorithm: 
1: if (St)x is not empty then 

St^(St\{(S t ) x }) 
end if 

S T <r- S T U {r = } 
return S 

Now we present the main algorithm. The general structure is as follows: In each iteration, a 
system S is selected from a list P of unfinished systems. An equation or inequation q is chosen from 
the queue Sq according to the selection strategy. Then q is reduced modulo St and incorporated 
into the candidate simple system St with the splitting algorithms as described above. In doing 
so, the algorithm may add new systems Si to P. As soon as the algorithm produces a system 
containing an equation c— for c G F \ {0} or the inequation 0^, this system is discarded. 

Algorithm 2.25 (Decompose). The algorithm is printed on page 13. 

We demonstrate the algorithm with a simple example. Note, that we will omit systems which 
are obviously inconsistent. 

,•2 



Example 2.26. Let S = (St,Sq) := (0,{(ar + x + 1)=, (x + a)^}) with a < x. According to 



Select, q := (x 2 + x + 1) = is chosen. As init(g) = 1 and reso(g, -S-q,x) = 1, the original system S 



is replaced by {{{x 2 + x + 1) = }, {(x + a)^}) . 

Now, q := (x + a)^ is selected and ResSplitDivide(5, (St)x, q) computes reso((ST)x,q,x) = 
prem((ST)x,q,x) = a 2 - a + I, resi((ST)x,q,x) = init(g) = 1, and res 2 ((ST)x,q,x) = 1. As St 
contains no equation of leader a, none of these polynomials can be reduced. Then, we decompose 
S into 

S := ( Kx 2 +z + l) = ,(a 2 -a + l)^M }j , 

=S T =Sq 
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Algorithm 2.25 (Decompose) 



Input: A system S' with (S') T = 0. 
Output: A THOMAS decomposition of S'. 
Algorithm: 

1: P<- {S"}; Result <- 

2: while |P| > do 

3: Choose SeF;F^P\{S} 

4: if |5 Q | = then 

5: Result «- EesuZi U {5} 

6: else 

7: g <- Select(S* Q ); 5 Q <- S Q \ {q} 
q <— Reduce(<?, St)', % ^- ld(q) 
if qi {0^, c = \ceF\ {0}} then 
if x 1 then 

if 5 is an equation then 

if (St)x is an equation then 

if Reduce(reso((S , T)x, q, %), St) = then 

(S,S u p) <- ResSplitGCD(5,q,a;); P^PU^} 
S <— lnsertEquation(S', p = ) 
else 

Sq<^ SqL> {q=, res ((S T ) x , q, %)=} 
end if 
else 

if (St)x is an inequation" then 

s Q ^s Q u {(s T ) x y, s t ^s t \ {(St)A 

end if 

{S, S 2 ) <- lnitSplit(5', g);Pf-PU {S 2 } 
(S,S 3 ,p) <- ResSplitSquareFree(S',g); P<-PU{5 3 } 
5 <— lnsertEquation(S l , p—) 
end if 

else if q is an inequation then 
if (St)x is an equation then 

(S,S 4 , P ) <- ResSplitDivide (S, (St)*,?); P^PU{S 4 } 
5 <— lnsertEquation(5,p = ) 
else 

(5, fife) <- lnitSplit(5, g);P^PU {S 5 } 
(S,S 6 ,p) <- ResSplitSquareFree(S*,g); P <- P U {S 6 } 
if (St)* is an inequation then 

(S,S 7 ,r) <- ResSplitDivide (S, (S T ) X , p); P <- PU{S T } 

(St)x <~ ( r 
else if (St)x is empty then 

(S T )x <~P=£ 

end if 
end if 
end if 
end if 

P^PU{S} 
end if 
end if 
end while 
return Result 



"Remember that (St)x might be empty, and thus neither an equation nor an inequation. 
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which is already simple, and 

Si := ({(.t 2 + x + 1) = }, {(x + a)?, (a 2 -a + 1) = }) 

=(Si)t =(Si)q 

We replace S\ by 



Si := ({{x 2 +x + 1) =) (a 2 - a + 1)=}, {(x + a) # }) 
de(Si,( 

/io/ds and Si is replaced with 



and apply ResSplitDivide(Si, ((Si)t)x> q) to Si again. This time, Reduce(a 2 — a + 1, (Si)t) = 



Si:=({ (x-a+l)= , (a 2 -a + !)=},{!#}) . 

pquo(a; 2 +x+l,x+a,cc) 

Finally, a THOMAS decomposition of S is: 

({{x 2 + x + 1) = , (a 2 - a + 1)^}, {(i - a + 1) = , (a 2 - a + 1)=}) . 

Proof (Correctness). First, note that it is easily verified that the input specifications of all 
subalgorithms are fulfilled (in particular, for lines 14 and 29, cf. Remark (2.7)(1)). 

The correctness of the Decompose algorithm is proved by verifying two loop invariants: 

1. P U Result is a disjoint decomposition of the input S' . 

2. For all systems S G P U Result, St is triangular and 

(a) 0<ir,a(p) is square-free and 

(b) </»a(init(p)) ± 

for all p e St with ld(p) = £ and all a e ©oI((St)<2; U (Sq) <x ). 

We begin with proving the first loop invariant. Assume that P U Result is a disjoint decom- 
position of S' at the beginning of the main loop. It suffices to show that all systems we add to P 
or Result add up to a disjoint decomposition of the system S, that is chosen in line 3. If Sq = 
holds in line 4, the algorithm just moves S from P to Result. 

In line 17, adding reso (( St ) a; , q, x) = to S does not change the solutions of S, as 4><x,a((ST)x) = 
and <jxx,a{q) = for each a 6 F[ y | y < x ] implies a (reso(p, q, x)) = (cf. (Mishra, 1993, 
Lemma 7.2.3)). 

Note now that if (S,Si) is the output of any of the ResSplitGcd, InitSplit, ResSplitSquareFree 
and ResSplitDivide algorithms, then (S U {q},Si) is a disjoint decomposition of So U {q}, where 
Sq is the input of the respective algorithm. It remains to be shown that the actions performed in 
lines 15, 25, 30, 36 and 38 are equivalent to putting q back into the system S. 

Let a G &ol(S <x ). In the context of line 15, Algorithm (2.19) guarantees 

0<x,a(p) = (f><x, a (( S T)x) = and (f> <x<a .(q) = . 
In the context of line 30, Algorithm (2.20) ensures that 

0<z,a(p) = 4=> (pKxAiS^x) = and (j) <x , a {q) ^ . 
In lines 25, 36 and 38, p has the same solutions as q, due to Algorithm (2.21) and 

/ x ^ ( t ) <x., a {q) <i><x, a {q) 

gcd(0 <2 , ia (g),(/> <:I ., a (^g)) gcd(0 <Xia (g), ^<?!> <X)a (g)) 

In addition, in line 36, 

<P<xA r ) ~ i / ■ %|ttt — nr =^ ^<*>a( r • p) ~ lcmO^aCCSr)*), <^ <!B ,a(») • 

gcd(0 <2;!a ((5 T ) :!; ), <a; ,a(p)) 
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This concludes the proof of the first loop invariant. 

Now, we prove the second loop invariant. At the beginning, the loop invariant holds because 
S' T = holds for the input system S'. Assume that the second loop invariant holds at the beginning 
of the main loop. 

One easily checks that all steps in the algorithm allow only one polynomial (St)x in St for 
each leader x, thus triangularity obviously holds. 

We show that all polynomials added to St have non-zero initial and are square-free. For 
©oI(S<a;) = 0, the statement is trivially true. So, let a £ &ol(S <x )- 

For the equation p— added as conditional gcd of (St)x an d q in line 15, it holds that 4><x,&(jp) 
is a divisor of 4><x,&{{St)x)- As 4><x,a((ST)x) is square-free by assumption, so is 4>< x ,a.(j>)- The 
inequation added to S in ResSplitGCD is by Definition (2.13) the initial of p—. 

The equation p— inserted into St in line 25 and the inequation p^ inserted in line 38 are 
square-free due to Algorithm (2.21) and their initials are non-zero as p is either identical to q, 
or it is a pseudo quotient of q by PRSi (g, -^q,x) for some i > 0. On the one hand, if p equals 
q, the call of InitSplit for q ensures a non-zero initial for p. On the other hand, the polynomial 
PRSi (g, -§^q, x) has initial res^ (g, J^q, x) , which is added as an inequation by ResSplitSquareFree. 
This implies that the initial of the pseudo-quotient is also non-zero. 

The equation p— that replaces the old equation (St)x in line 30 is the quotient of (St)x by an in- 
equation. It is square-free, because 4>< x ,bSj>) is a divisor of <fi<x,a((ST)x), which is square-free by as- 
sumption. Again, p is either identical to (St)x or a pseudo quotient of (St) x by PRSi ((St)x, q, x) 
for some i > and, using the same arguments as in the last paragraph, the initial of p does not 
vanish. 

Finally, consider the inequation (r-p)^ added in line 36 as a least common multiple of ((St)x)^ 
and p^. The inequation </><£, a(p) is square-free and has non- vanishing initial for the same reasons 

as before. Due to <j> <Xt *(r) ~ gcd^J^Csr^U^Ap)) ' the P oh / nomials 4><x,a( r ) and < P<yAp) 
have no common divisors. As <j)<cx,a{i') divides 4>< x ^{{St) x ), using the same arguments as before, 
4'<x,a(j') is square-free and has a non-vanishing initial. This completes the proof of the second 
loop invariant. 

It is obvious that a system S with Sq = for which these loop invariants hold is simple. Thus 
the algorithm returns the correct result if it terminates. □ 

We now start showing termination. The system S chosen from P is treated in one of three 
ways: It is either discarded, added to Result, or replaced in P by at least one new system. To 
show that P is empty after finitely many iterations, we define an order on the systems and show 
that it is well-founded. Afterwards we prove termination by detailing that the algorithm produces 
descending chains of systems. 

Definition 2.27. For transitive and asymmetric 1 partial orders <j for i = l,...,m, we define 
the composite order " < " := [<i, . . . , < m ] as follows: a < b if and only if there exists i G 
{l,...,m} such that a <i b and neither a <j b nor b <j a for j < i. The composite order 
is clearly transitive and asymmetric. An order < is called well-founded, if each <-descending 
chain becomes stationary. 

The following trivial statement will be used repeatedly: 

Remark 2.28. If each < t is well-founded, then so is the composite ordering <, using the notation 
from Definition (2.27). 

Now we define the orders and show their well-foundedness: 

Definition and Remark 2.29. Define -< as the composite order -<2, -<3, -<&] of the four or- 
ders defined below. It is well-founded since the -<j are. 



7 A relation H is asymmetric, if S -< 5" implies 5" S for all S,S'. Asymmetry implies irreflexivity. 
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1. For i — l,...,n define <\. Xi by S <i. Xi S' if and only if mdeg ((St)x) < mdeg ((S T ) Xi ) , 
with mdeg ((St)s-) '■— 00 */ (St)^- is empty. Define the composite order -<i as [-<i, Kl 
, • • • j ~<i,x„]' Since degrees can only decrease finitely many times, the orders -<\ >Xi are clearly 
well-founded and, thus, -<i is. 

2. Define the map fi from the set of all systems over R to {1, X\, . . . , X n , x^}, where fJ,(S) 
is minimal such that there exists an equation p 6 (Sq)^?^ with Reduce(S l T,p) 7^ 0, or 
n(S) = Xoo if no such equation exists. Then, S S' if and only if fJ-(S) < fJ-(S') with 
1 < Xi and Xi < Xoq for i E {1, . . . ,n}. The ordering -<2 is well-founded since < is well- 
founded on the finite set {1, xi, . . . , x n , Xoo}. 

3. S -<3 S' if and only if there is p-t E R* and a finite (possibly empty) set L C R* with 
ld(g) < ld(p) V q G L such that Sq l±) {p^} = S'q W L holds. We show well-foundedness by 
induction on the highest appearing leader x in (Sq)*: Forx = l we can only make a system 
S -<3~smaller by removing one of the finitely many inequations in (Sq)*. Now assume that 
the statement is true for all indeterminates y < x. By the induction hypothesis we can 
only -<3~decrease S finitely many times without changing (Sq)* . To further ^-decrease S, 
we have to remove an inequation from (Sq)* . As (Sq)* is finite, this process can only be 
repeated finitely many times until (Sq)* = 0. Now, the highest appearing leader in (Sq)* is 
smaller than x and by the induction hypothesis, the statement is proved. 

4- S ~<i S' if and only if \Sq\ < \S'q\. 

Proof (Termination). We will tacitly use the fact that reduction never makes polynomials 
bigger in the sense of Remark (2.7) (3). 

We denote the system chosen from P in line 3 by S and the system added to P in line 43 by S. 
We prove that the systems S, S%, . . . , S7 generated from S are ^-smaller than S. For i — 1, ... ,4 
we will use the notation S ^ . S' if neither S <i S' nor S' <i S holds. 

For j — 1, . . . , 7, ((Sj)r) = — (St) = and thus Sj ^ S. The properties of Select in Definition 

(2.22) directly require, that there is no equation in (Sq) = with a leader smaller than x. However, 
the equation added to the system Sj returned from InitSplit (2.12) is the initial of q, which has a 
leader smaller than x and does not reduce to (cf. Remark (2.7)(2)). Furthermore, the equations 
added in one of the subalgorithms based on ResSplit (2.18) have a leader smaller than x and do 
not reduce to 0. In each case Sj -<2 S is proved. 

It remains to show S -< S. If g is reduced to = , then it is omitted from Sq and so S -<4 S. As 
the system is otherwise unchanged, S ^ . S, i — 1,2,3 and therefore S -< S holds. If q is reduced to 

for some c £ F \ {0}, then S -<3 S and S B. S,i = 1, 2, since the only change was the removal 
of an inequation from Sq. Otherwise, exactly one of the following cases will occur: 

Lines 14-15 set (St)x to p = of smaller degree than (St)x and 20-25 add (St)x as a new 
equation. In both cases we get S -<i S. 

In line 17, St = St implies S ^ S. The polynomial q is chosen according to Select 

(cf. (2.22)(1)), which implies (Sq)< x = and (Sq)< x = {reso((ST) x , q, %)=}■ Line 13 ensures 
Reduce(reso((Sr)a:, q, x), S) ^ and, thus, S -<2 S follows. 

Consider lines 29-30. If the degree of (St)x is smaller than the degree of (St)x, then S -<% S. 
In case the degree doesn't change, we have S ^ S and (Sq) = = (Sq) = guarantees S ^ 2 S. 
However, q is removed from Sq and replaced by an inequation of smaller leader, which implies 
S^ 3 5. 

In 31-39, obviously S ^. S, i = 1,2. As before, q is removed from Sq and replaced by an 
inequation of smaller leader, which once more implies S -<3 S. □ 
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2.3. Notes to Applications of Simple Systems 

In this subsection, we shortly present some examples where simple systems are necessary and 
any weaker decomposition into triangular systems is not sufficient. 

The properties of simple systems (cf. Definition (2.2)) correspond exactly to the following 
fibration structure on the solution sets (cf. Plesken (2009a)). Let S be a simple system and 
Hi : F l —} F l : (en, . . . ,a;) i-» (oi, . . . , a$-i)- Furthermore, for any solution a £ &ol(S< Xi ), let 
S , a = nr 1 ({n,(a)}). Then, if S Xi is an equation, |sj a | = nl dcg(S , 2;i ) holds. If S Xi is an inequation, 
then Si a — F \ s^ a with |Sj. a | = mdeg(5 :Ei ). If S Xi is empty, then s i a = F. The cardinalities of 
s«, a or Sj ja are constant for each i, i.e. independent of the choice of the solution a £ 6o[(S< Xj ) (cf. 
Remark (2.3)). We can examine solution sets of arbitrary systems by decomposing them disjointly 
into simple systems. Further analysis of this fibration structure, especially in the context of 
algebraic varieties, is a topic of future research. 

We already saw such a fibration structure in Example (2.1). In this case, other triangular 
decompositons like a decomposition into regular chains would have only resulted in a single system 
consisting of the polynomial p from the input. 

A special case occurs when all polynomials in the input and output can be factored into linear 
polynomials. If we compute the counting polynomial as introduced by Plesken (2009a) (which 
requires the disjointness of the decomposition and the fibration structure), we can substitute the 
cardinality of a finite field F (of sufficiently large characteristic) into the counting polynomial of a 
Thomas decomposition computed over Q. This yields the exact number of distinct solutions over 
F. For example, the counting polynomial of a THOMAS decomposition of {det(M)^} for a generic 
n x n matrix M = (xij)i<<j<n yields the well-known formula for the cardinality of GL n (F) for 
any finite field F. Furthermore, we can automatically reproduce the results in (Plesken, 1982, 
Ex. V.4), where pairs of matrices (A, B) with given ranks of A, B, and A + B are counted. 

Plesken (2009b) gave another example concerning the GAUSS-BRUHAT-decomposition and the 
LU-decomposition. The cells of these decompositions of M as above can be identified with certain 
simple systems in the THOMAS decomposition of {det(M)^} for suitable rankings on the Xij. 

We clearly see that simple systems are necessary for these applications to expose the aforemen- 
tioned fibration structure and count solutions. A disjoint decomposition into triangular systems 
with weaker properties does not suffice. 

3. Differential Thomas Decomposition 

The differential Thomas decomposition is concerned with manipulations of polynomial dif- 
ferential equations and inequations. The basic idea for our construction of this decomposition is 
twofold. On the one hand, a combinatorial calculus developed by Janet finds unique reductors 
and all integrability conditions by completing systems to involution. On the other hand, the al- 
gebraic Thomas decomposition makes the necessary splits for regularity of initials and ensures 
disjointness of the solution sets. 

Initially, we recall some basic definitions from differential algebra. Then, we summarize the 
Janet division and its relevance. Its combinatorics lead us to substitute the algebraic algorithm 
InsertEquation by its differential analog. Afterwards, we review a differential generalization of the 
algebraic reduction algorithm and present the algorithm Reduce utilized for differential reduction. 
Replacing the insertion and reduction from the previous section with these differential counterparts 
yields the differential THOMAS decomposition algorithm. 

3.1. Preliminaries from Differential Algebra 

Let A = {d\, . . . , d n } be a non-empty set of derivations and F be a A-ring. This means any 
dj £ A is a linear operator dj : F —> F which satisfies the Leibniz rule. Given a differential 
indeterminate u, the polynomial A-ring F{u} := F [ u\ | i € Z™ ] is defined as the polyno- 
mial ring infinitely generated by the algebraically independent set (m)a '•= {Ui | i £ Z" }. The 
operation of dj £ A on (it) a is defined by djUi = u- 1+ej and this operation extends linearly and 
via the Leibniz rule to F{u}. Let U = {u^, . . . ,M^ m ^} be a set of differential indeterminates. 
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The multivariate polynomial A-ring is given by F{U} := FjV 1 '} . . . {u< m )}. Its generators, the 
elements of (U)a '■= | i G €{!>••• 7 TO }|i are called differential variables. From 

now on let F be a computable A-field of characteristic zero. 

The differential structure of F uniquely extends to the differential structure of its algebraic 
closure F (Kolchin, 1973, §11.2, Lemma 1). Let E := (B™ =1 F[[ * • * j Zn}] where ^[f^i, . . . , z n ]] 

denotes the ring of formal power series in zi, ... , z n . Then E is isomorphic to F ^ A via 
a:®F[[ Zl ,...,z n ]} ->¥ 

3=1 

where z 1 :— z 1 ^ ■ . . . ■ z l ™ and i! := i%\ ■ ■ ■ ■ • 

We define solutions in E, consistent with the algebraic case: For e G E, let 

4> e : F{U} -> F : up h> a(e)(tti (i) ) 

be the F-algebra homomorphism evaluating the differential variables at e. A differential equa- 
tion or inequation for m functions U — {u^\ . . . , i*( m )} in n indeterminates is an element 
p G F{U}, written p = or p^, respectively. A solution of p = or p^ is an e G E with <p e ((p)A) = {0} 
or (f>e({p)A) 7^ {0}, respectively. Here, (p)& denotes the differential ideal in F{U} generated by p. 
Furthermore, e G E is called a solution of a set P of equations and inequations, if it is a solution 
of each element in P. The set of solutions of P is denoted by &ol(P) :— &oIe{P) Q E. 

In differential algebra one usually considers solutions in a universal A-field, while we consider 
power series solutions. As the universal differential field we can take the universal closure F of F. 
There is a strong link between these two concepts. On the one hand, Seidenberg (1958, 1969) has 
shown that every finitely differentially generated differential field is differentially isomorphic to a 
differential field of meromorphic functions in n variables. On the other hand, . . . , z n ]} ^ 

F ((zi, . . . , z n )) <— > F. Here, the first map is the natural embedding into the quotient field and the 
second is an embedding given by the definition of the universal A-field (Kolchin, 1973, §11.2 and 
§111.7), as F((zi, . . . , z n )) is a finitely generated A-field extension of F. Thus, any power series 
solution can be considered as a solution in the universal differential field. 

A finite set of equations and inequations is called a (differential) system over F{U}. We will 
be using the same notation for systems as in the algebraic Thomas decomposition introduced in 
§2.1 and §2.2, in particular a system S is represented by a pair (St, Sq). However, the candidate 
simple system St will also reflect a differential structure based on the combinatorics from the 
following section. 

3.2. Janet Division 

In this subsection we will focus on a combinatorial approach called Janet division (cf. Gerdt and Blinkov 
(1998a)). It manages the infinite set of differential variables and guarantees inclusion of all inte- 
grability conditions in a differential system. For this purpose, it partitions the set of differential 
variables into "free" variables and finitely many "cones" of dependent variables. We present an 
algorithm for inserting new equations into an existing set of equations and adjusting this cone 
decomposition accordingly. An overview of modern development on Janet division can be found 
in Gerdt (2005) and Seiler (2010) and the original ideas were formulated by Janet (1929). 

A (differential) ranking < is defined as a total order on the differential variables and 1 with 
1 < u V u G U, such that 

1. u < djU and 

2. u < v implies djU < djV 

for all u, v G (U)a, dj £ A. From now on let < be an arbitrary and fixed differential ranking. For 
any finite set of differential variables, a differential ranking induces a ranking as defined for the 
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algebraic case in §2.1. Thereby, in accordance to the algebraic part, define the largest differential 
variable \d(p) appearing in a differential polynomial p £ F{U} as leader, which is set to 1 for 
p £ F. Furthermore, define mdeg(p) and init(p) as the degree in the leader and the coefficient of 
ld(p) mdog ( p \ respectively. 



Example 3.1 



Consider two derivations A = {d x ,d t } and one differential indeterminate u. 

In this setting, any partial differential equation with constant 
coefficients in one dependent variable and two independent vari- 
ables can be represented as a differential polynomial in C{u}. 

The ranking < is defined by Ui lt i 2 < Uj ± j 2 if and only if either 
i\ + t2 < ji + h or ii + i 2 = ji + ]2 and i 2 < ji holds. Thus, 
the smallest differential variables are: u , < u 1-0 < u 01 < u 2Q < 
1*1,1 < u 02 < u 3 a . When considering the set of differential vari- 
ables as a grid in the first quadrant of a plane, the picture on the 
left illustrates this ranking. 

Consider (u OA +u 0:0 u 1} o)= representing the inviscid Burger's 
0. As in the algebraic part, we indicate an equation in the picture by 
attaching it to its leader. However, contrary to the algebraic part, a differential equation does 
not only affect its leader, but also the derivatives of its leader. This is because property 2 of a 
differential ranking implies (91d(p) = \d(dp) V<9 £ A,p6 F{U}. For example dt{u o l + u 00 u l a ) = 
Ua,2 + u 01 u l a + u o a u 11 . In the diagram we illustrate this by drawing a cone with apex u 01 . 




equation ^ + u|| 



(U ,i + Mo.oU li0 ) 



(u 2 




Assume that we are only interested in solu- 
tions of the inviscid BURGER 's equation which 
are linear in x. So, we add the second equa- 
tion (it 2 ,o)= to our system. This second equa- 
tion also affects the derivatives of its leader. In 
particular, {u QA +u QS ,u l a ) = and (u 2 ) = both af- 
fect the differential variable u 2 l and its deriva- 
tives. This contradicts the triangularity of the 
system. According to the involutive approach 
as suggested by Janet, we don't allow certain 
equations to be derived by certain partial derivations. In this example, we allow (u 2 ,o)= to be de- 
rived only by d x . In the diagram we illustrate this by drawing a (degenerate) cone with apex u 2 Q 
in direction of d x . Thus, the differential consequence (<9tU 2 ,o)= is not yet considered and, so, we 
have to add it as a separate equation for further treatment. 

A set W of differential variables is closed under the action of A' C A if dtw £ W for all 
di £ A' and w £ W. The smallest set containing a differential variable w, which is closed under 
A', is called a cone and denoted by {w)a'- In this case, we call the elements of A' reductive 
derivations 8 . The A'-closed set generated by a set W of differential variables is defined as 



(w) A > ■■= n w i ^ v)* 



W i DW 

Wi A'-closed 



For a finite set W — {wi, . . . ,w r }, the Janet division algorithmically assigns reductive 
derivations to the elements of W such that the cones generated by the w £ W are disjoint 
(cf. Gerdt et al. (2001) for a fast algorithm). We call these derivations jANET-reductive. The 
derivation di £ A is assigned to the cone generated by w = iij € W as reductive derivation, if 
and only if 

i ; = max ji; | u[P £W,'i' k — ife for all 1 < k < Z j 



8 In Gerdt (1999) and (Seller, 2010, Chap. 7) the reductive derivations are called multiplicative variables and in 
Bachler et al. (2010) they are called admissible derivations. 
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holds (Gerdt, 2005, Ex. 3.1). We remark that j is fixed in this definition, i.e., when constructing 
cones we only take into account other differential variables belonging to the same differential 
indeterminate. Furthermore, the assignment of reductive derivations to w £ W in general depends 
on the whole set W. The reductive derivations assigned to w are denoted by Aw{w) C A and 
we call the cone (w)a w ( w ) the Janet cone of w with respect to W. This construction ensures 
disjointness of cones but not necessarily that the union of cones equals (W)a- The problem 
is circumvented by enriching W to its Janet completion W 3 W. This completion W is 
successively created by adding any 

w = d l w j <£ |+J {w) Aw{w) 

to W, where Wj £ W and 9j£A \ A^(wj). This leads to the disjoint Janet decomposition 

(W) A = Wj™>A wW 

that algorithmically separates a A-closed set (W)a into finitely many cones (w)a~ (to)- For details 
see (Gerdt, 2005, Def. 3.4) and (Gerdt and Blinkov, 1998a, Cor. 4.11). 

We extend the Janet decomposition from differential variables to differential polynomials ac- 
cording to their leaders. To be precise, At(<?) := A ld( - T )(ld(q)) for finite T C F{U} and q £ T. 
We call a derivative of an equation by a finite (possibly empty) sequence of derivations a prolon- 
gation. If all these derivations are reductive, the derivative is called reductive prolongation of 
q with respect to T. Otherwise it is called non-reductive prolongation. 

A differential polynomial p £ F{U} is called reducible modulo q £ F{U}, if there exists 
i G Z£ such that d[ 1 -...-d i r ?\d(q) = ld^ 1 . .-d%q) = ld(p) and mdeg(^ 1 • . . . • d% q) < mdeg(p). 
For i 7^ (0, . . . , 0) the condition on the main degree always holds. We now restrict ourselves to 
reductive prolongations: For a finite set T C F{U}, we call a differential polynomial p £ F{U} 
jANET-reducible modulo q £ T w.r.t. T, if p is reducible modulo q and c^ 1 • . . . • d^q is a reductive 
prolongation of q w.r.t. T, with i £ Z™ from the reducibility conditions. We also say that p is 
JANET-reducible modulo T if there is a q £ T such that p is JANET-reducible modulo q w.r.t. T. 

A set of differential variables T c (U)a is called minimal, if for any set S C (U)a with 

A T (t) = |+J(s)a s ( s ) 

tGT s£S 

the condition T C S holds (Gerdt and Blinkov, 1998b, Def. 4.2). We call a set of differential 
polynomials minimal, if the corresponding set of leaders is minimal. 

At each step of the algorithm we assign reductive derivations to the equations in (St) = - 
When an equation p is not reducible modulo (St) = , it is added to (SV) = - Then, we remove 
all polynomials from St that have a leader which is derivative of ld(p). This will later ensure 
minimality. In addition, when adding a new equation to (St) = , all non-reductive prolongations 
are put into the queue. This is formalized in the following algorithm. 

Algorithm 3.2 (InsertEquation). 

Input: A system S' and a polynomial p = £ F{U} not reducible modulo (S' T ) = . 
Output: A system S , where (St) = Q {S't) = U {p=} * s maximal satisfying 

(ld(Sr) \ {ld(p)» n (ld(p)) A = 0, 
S Q = S' Q U (S' T \ S T ) U {{d iq ) = | q £ (SrTA ? A ((Sr) =)( 9 )} ■ 

Algorithm: 
1: S^- S' 
2: St <— St U {p=} 
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3: for q G St\ {p} do 

4: i/ld(g) G (ld(p)) A then 

5: Sq^SqU {q} 

6: S T <- St\ {q} 
7: end if 

8: end for 

9: Reassign reductive derivations to (St) = 

10: S Q ^-S Q U {(CW)= I q G (5^)=,$ £ A ((St) =) (<?)} 

11: return S 

Correctness and termination are obvious. We remark that a non-reductive prolongation might 
be added to Sq several times. An implementation should remember which prolongations have 
been added before to avoid redundant computations. 

3.3. Differential Simple Systems 

In this subsection, we extend the algebraic reduction algorithm to its differential counterpart. 
Finally, we can define differential simple systems at the end of this subsection. 

The Janet partition of the dependent differential variables into cones provides a mechanism 
to find the unique reductor for the differential reduction in a fast way (cf. Gerdt et al. (2001)). 
We prolong this reductor and afterwards apply a pseudo reduction algorithm. 

For a valid pseudo-reduction, we need to ensure that initials (and initials of the prolongations) 
of equations are non-zero. Let r G F{U} with x = ld(r) and define the separant sep(r) := 
One easily checks that the initial of any non-trivial prolongation of r is sep(r) and the separant 
of any square-free equation r is non-zero (cf. (Kolchin, 1973, §1.8, Lemma 5) or (Hubert, 2003b, 
§3.1)). So, by making sure that the equations have non- vanishing initials and are square-free, as 
in the algebraic case, we ensure that we can reduce modulo all prolongations of r. This provides 
the correctness of the following reduction algorithm. 9 

Algorithm 3.3 (Reduce). 

Input: A differential system S and a polynomial p G F{U}. 

Output: A polynomial q that is not J ANET -reducible modulo St with <f> e {p) = if and only if 

(j> e (q) = for each e G &ol(S). 

Algorithm: 

1: X <r- ld(p) 

2: while exists q = G (5*t) = and i G Z" with ij = for dj ^ A(s r )=(q) such that • . . . • 

<9 r \" ld(q) = ld(p) and mdeg(9} 1 • . . . • d]?p) > mdeg(q) hold do 
3: p <- prem (p, d} 1 • . . . ■ d% q, x) 
4: x <— ld(p) 
5: end while 

6: if Reduce^, init(p)) = then 

7: return Reduce(S> - init(p)a; mdcg W) 

8: else 

9: return p 
10: end if 

A polynomial p G F{U} reduces to q modulo St if Reduce(5,p) = q. A polynomial 
p G F{U} is called reduced 10 modulo St if it reduces to itself. The properties of the algebraic 
reduction algorithm from Remark (2.7) also apply for this reduction algorithm. 



9 In differential algebra, one usually distinguishes a (full) differential reduction as used here and a partial 
(differential) reduction. Partial reduction only employs proper derivations of equations for reduction (cf. (Kolchin, 
1973, §1.9) or (Hubert, 2003b, §3.2)). This is useful for separation of differential and algebraic parts of the algorithm 
and for the use of Rosenfeld's Lemma (cf. Rosenfeld (1959)), which is the theoretical basis for the Rosenfeld- 
Ghobner algorithm (cf. Boulier et al. (2009, 1995); Hubert (2003b).) 

10 There is a fine difference between not being reducible and being reduced. In the case of not being reducible 
the initial of a polynomial can still reduce to zero and iteratively the entire polynomial. 



21 



Termination of the reduction algorithm is provided by Dickson's Lemma (cf. (Cox et al., 
1992, Chap. 2, Thm. 5) or (Kolchin, 1973, §0.17, Lemma 15)), which states that the ranking < is 
well-founded on the set of leaders, i.e., a strictly <-descending chain of leaders is finite. 

Example 3.4. We continue Example (3.1) and take care of the differential consequence (u 2 ,i)=- 
We reduce (it 2 ,i)= modulo the system S with 



St := \pi ■= (mq,i + Hq,oMi,q)=,P2 : = (u 2 , )=\ 




odulo p 2 modulo d x P2 



First, we observe, that ld(u 2 ,i) = u 2 ,i is in the cone gen- 
erated by ld(pi) and \d(d^pi) = ld(u 21 ). Thus, we re- 
duce (u 21 ) modulo d^Pi and the pseudo reduction yields 
« 0i0 ti 3i0 + 3u li0 u 2>0 . Second, we reduce u„ u 3 , + 3u li0 ti 2 ,o 
modulo d x p2, because ii 3)0 lies in the cone generated by 
{u2,o)=- This results in 3u lt0 u 2 , and a third reduction step 
modulo p2 produces zero. As a result, the only differential consequence is already implied by the 
system. In this desirable situation, there are no further integrability conditions, which motivates 
the definition of involutivity below. 

Now, we define differential simple systems. We demand algebraic simplicity, involutivity of 
differential equations as seen in the previous Example (3.4), and some minimality conditions. 

Definition 3.5 (Differential Simple Systems). A differential system S is ( Janet-,) involu- 
tive, if all non-reductive prolongations of (St) = reduce to zero modulo (St) = ■ 
A system S is called differentially simple or simple, if 

1. S is algebraically simple (in the finitely many differential variables that appear in it), 

2. S is involutive, 

3. S = is minimal, 

4- no inequation in S^ is reducible modulo S = . 

A disjoint decomposition of a system into differentially simple subsystems is called (differential) 
Thomas decomposition. 

As in the algebraic case, every simple system has a solution in E. 

3.4- Differential Decomposition Algorithm 

The differential Thomas decomposition algorithm is a modification of the algebraic Thomas 
decomposition algorithm. We have already introduced the new algorithms InsertEquation (3.2) 
for adding new equations into the systems and Reduce (3.3) for reduction, that can replace their 
counterparts in the algebraic algorithm. This subsection provides the necessary correctness and 
termination proofs for this modified algorithm. It then demonstrates this algorithm with examples. 

Algorithm 3.6 (DifferentialDecompose). 

Input: A differential system S' with (S')t = 0- 
Output: A differential THOMAS decomposition of S' . 

Algorithm: The algorithm is obtained by replacing the two subalgorithms InsertEquation and Reduce 
in Algorithm (2.25) with their differential counterparts (3.2) and (3.3), respectively. 

Proof (Correctness). The correctness proof of the algebraic decomposition Algorithm (2.25) 
also holds verbatim for the differential case. Therefore, we do not need to show that the output 
is algebraically simple. We will prove three loop invariants for any system S £ P U Result: 
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1. (<St) = is minimal. 

2. No inequation in (St)^ is JANET-reducible modulo St- 

3. Let r be any non-reductive prolongation of (St) = ■ Then r reduces to zero by using both con- 
ventional differential reductions 11 of (Sq) = and reductions modulo reductive prolongations 
of (StT- 

The first loop invariant is a purely combinatorial matter, which is proved by Gerdt (2002) for 
an algorithm using exactly the same combinatorial approach. 

Proving the second loop invariant is equally simple. On the one hand, a newly added inequation 
q in St is not JANET-reducible modulo (St) = , since algorithm Reduce (3.3) is applied to it before 
insertion. On the other hand, algorithm InsertEquation (3.2) removes all inequations from St 
which are divisible by a newly added equation and places them into Sq. 

The third loop invariant clearly holds at the beginning of the algorithm, because St is empty. 

We claim that reduction of an equation q— £ Sq by (St) = in line 8 of Algorithm (2.25) does 
not affect the loop invariant, i.e. any non-reductive prolongation r reducing to zero beforehand 
reduces to zero afterwards. We prove this claim by performing a single reduction step on q, which 
generalizes by an easy induction. Let q' := prem(q,p, x) = to • q — pquo(q,p, x) ■ p be a pseudo 
remainder identity (see (1) on page 5) reducing q to q' modulo p. Then a pseudo remainder 
identity prem(r, q, x) = m' ■ r — pquo(r, q,x) ■ q describing a reduction of r modulo q might simply 
be rewritten as the iterated identity 

m ■ prem(r, q, x) — to ■ m ■ r — pquo(r, q, x) ■ q' — pquo(r, q, x) ■ pquo(g,p, x) ■ p. 

prem(prem(r,p,a:),Q / ,x) 

Using the Leibniz rule the same holds for reduction modulo partial derivatives of q. This holds 
especially for an equation q— £ Sq reducing to modulo (5*t) = in line 8, which can be removed 
from Sq without violating the loop invariant. 

Now, we consider line 25, where InsertEquation inserts the square-free part p = of q = into St and 
show that this does not violate the third loop invariant. First, the non- reductive prolongations in 
{(<?,r) = | r £ [ST) = di 4- ^((St)=) ( r )l are added to Sq as equations. Thus, any of these reduce to 
modulo (Sq) = . Second, moving equations from St back into Sq in InsertEquation does not change 
the loop invariant either, because their reductive prolongations can still be used for reduction 
afterwards. Third, every non-reductive prolongation that reduced to zero using q = £ (Sq) = still 
reduces to zero after InsertEquation. This holds for two reasons. On the one hand, everything that 
reduces to zero modulo q = , also reduces to zero modulo p = . Write to • q = p ■ qi with ld(m) < x 
and 4> a (to) ^ Va £ ©o[(S' < id(g))- Then p algebraically pseudo-reduces q to zero. Any derivative 
dq of q is reduced to zero modulo p— and (<9p) = , since d(m ■ q) = (dp) ■ q\ +p - (dqi) for any d £ A. 
Inductively, the same holds for repeated derivatives of g = . Therefore, p— implies all constraints 
given by q=. On the other hand, all reduction steps modulo p— are either Janet-reductions modulo 
p— w.r.t. St or differential reductions modulo non-reductive prolongations of p = . The latter 
equations have been added to Sq. 

When computing the gcd of two equations in line 14, the gcd of q and (St)% will be inserted 
into St and reduces everything to zero that both q and (St)x did. As above, the non-reductive 
prolongations are covered by inserting them into Sq and the reductive prolongations are implied. 

Dividing an equation (St)x by an inequation in lines 29 and 30 also influences (St) = - 
The new equation p =1 being a divisor of (St) Xi reduces everything to zero that (St)x and its 
non-reductive prolongations did by the same arguments as before. 

This proves the third loop invariant. When the algorithm terminates, Sq is empty and thus 
all non-reductive prolongations from (5*t) = JANET-reduce to zero modulo (St) = - The system is 
therefore involutive. 



11 i.e. modulo any prolongation 
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Furthermore, the first loop invariant implies minimality and the second loop invariant implies 
that no inequation is reducible by an equation, since for an involutive set reducibility is equivalent 
to jANET-reducibility. □ 



Our main tool for proving the termination of the algorithm is using six orders on differential 
systems. These are similar to the four orders used to show the termination of the algebraic 
decomposition algorithm. We use Dickson's lemma as main tool to show the well-foundedness of 
these orders. 

Definition and Remark 3.7. Define the orders -<i a , -<\b, <ic, <2, -<3, and -<4 as follows. 

-<l a : For V C {U)a there is a unique minimal set v(V) C V with V C {v{V))a ( Cox et ai, 1992, 
Chap. 2, §4, exercise 7 and 8), called canonical differential generators of V . For a 
system S, define v(S) as i/(ld((Sr) = ))- For systems S,S' we define S -^\ a S' if and only if 
min < (i/(S')\i/(S")) < min < (z/(S")\^ , (5')). An empty set is assumed to have Xoo as minimum, 
which is <-larger than all differential variables. By Dickson's lemma, -<i a is well-founded. 

-<lb: For systems S, S' define S -<u S' if and only if S ^ S' and min< (ld((Sr) = ) \ ld((S"T) = )) < 
min< (ld((S"T) = ) \ ld((Sr) = )). Minimality of (St) = at each step of the algorithm and the 
constructivity property of the JANET division (Gerdt and Blinkov, 1998a, Prop. 4-.1S) imply 
well-foundedness o/-<i6 (Gerdt and Blinkov, 1998a, Thm. 4-H). 

<\ c : For systems S and S' with S ^ S' and S ^ S' , both (St) = and (S't) = have the same 
leaders X\,...,Xi. Define S -<i c ,x k S' if and only if mdegtXSV)^) < radeg((S T )^.). This 
order is clearly well-founded. For these systems define S -<\ c S' as [-<i ca:i , . . . , -Kicxjj which 
is again well-founded as a composite order. 

-<2-' This is defined identical to the algebraic -<2- We remark, that in this case the set of possible 
leaders is {1} U (U)a- To show well-foundedness of the differential ordering -<2 we use that 
< is well-founded on the set of leaders as implied by DlCKSON 's lemma. This way, < is 
extended to a well-founded ordering on {l^x^} U (£/)a with 1 < y and y < x^ for all 

ye(U) A . 

-< 3 ; This is verbatim the same condition and proof of well-foundedness as in the algebraic case. 
However, in the latter proof we do a NOETHERian induction (Bourbaki, 1968, III. 6. 5, 
Prop. 7) instead of an ordinary induction. 

-<4/ This is identical to the algebraic case. 

Remark (2.28) provides the well-foundedness of the composite order -<:= [^i a ,^ib, ^10^2, ^3,^4] • 

Proof (Termination). We prove termination the same way as in the algebraic case. All argu- 
ments where systems get -<2, ~<3, or -(4 smaller apply verbatim here. 

In the algebraic case a system -<i-decreases if and only if either an equation is added to St or 
the degree of an existing equation in St is decreased. We adapt this argument to the differential 
case: On the one hand, inserting a new equation with a leader that is not yet present in ld((5T) = ) 
decreases either -K\ a or -<ib. On the other hand, if an existing equation in (<S*t) = is replaced by 
one with the same leader and lower degree, the system -<i c -decreases. 

Thus, like in the algebraic termination proof, we have a strictly decreasing chain of systems 
and, thus, termination is proved. □ 

In the following examples, we use jet notation for differential polynomials, e.g., u XtXtV := U2.1 
in the case A = {d x ,d y } and U — {u}. 

We give an example taken from (Buium and Cassidy, 1999, pp. 597-600): 
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Example 3.8 (Cole-Hopf Transformation). For F := R(x, t), A = { JL, jfc}, and U = {77, C} 
consider the heat equation h = (rj t + T] xx )— and Burger's equation 6 = (Ct + Cxx + %Cx 1 C)=- 

First, we claim that any power series solution for the heat equation with a non-zero constant 
term can be transformed to a solution of BURGER 's equation using the Cole-Hopf transformation 
A : rj h4 — . A differential THOMAS decomposition for an orderly ranking with £ x > r\ t of 

{h = , (77 • C - »fe)=,»7/} 

V » ' 

consists of the single system 

S= {(Vx- v ()=, (v ■ Cx+ Vt + v C 2 )=,v_ i } 

and one checks that Reduce(5 f , b) = holds. This implies that A maps any non-zero solution of the 
heat equation to a solution of Burger 's equation. 

In addition we claim that A is surjective. For the proof we choose an elimination ranking ( cf. 
(Hubert, 2003b, §8.1) or Boulier (2007)) with 77 > (, i.e., 771 > <j for all i,j G Z> . We compute 
a differential THOMAS decomposition of {h—,b—, (17 • C — r) x )=,r)jt} . It consists of the single system 

S={(r h -v()=,(v(x+Ih+V( 2 U,b=X^} . 

The elimination ordering guarantees that the only constraint for £ is Burger's equation b—. As S 
is simple, for any solution f £ ©ol(&=) there exists a solution (g,f) G 6o[(5 l ) (cf. (2.3)), implying 
that A is surjective. 

Elements of the A-field F are not subjected to splittings and assumed to be non-zero. However, 
we are able to model the elements of F as differential indeterminates. For example for F = C(x) 
with A = {^}, we can study a differential polynomial ring over C{X} instead and replace x by X 
in all equations and inequations. We subject X to the relation -§^X — 1 for X being "generic" or 
{-§^X — 1) • -^X = 0, if we allow specialization of X. Both these cases are considered in examples 
(3.9) and (3.10), respectively, and will be subject of further study. 

Example 3.9. For F := C(x), A = Jj} and U — {u} consider the special case 

(iH - u xx - x ■ u x - u) = (5) 

of the Fokker-Planck equation. We add an auxiliary differential indeterminate X to U and 
instead examine the equation 

(u t ~ u xx ~ X ■ u x - u) = , (X x - 1) = , (X t ) = (6) 

in the A-ring <C{X,u}. An elimination ranking X 3> u splits the system (6) into two simple 
systems: 

(i) (u x ■ {~ u xxx + u xt - 2u x ) - u xx ■ (u t - u xx - u)) = , 
{u x ■ {- u xxt + utt - u t ) - u xt ■ (u t - u xx —u)) = , 
(u t - u xx - X ■ u x - u) = , (y ^ 

(ii) (u x ) = , (ut ~ u) = , (X* - 1) = , (Xt) = 

Due to the ranking, the first two equations in (i) generate (F{u}[A\ ■ (u t — u xx — X ■ u x — u)) D 
C{u}, i.e., they have constant coefficients. These two equations are the derivatives of M *~^^~ M _ Xj 
which is clearly equivalent to (5) in the case u x 7^ 0. 

The next example sketches an approach to treat equations with variable coefficients and find 
submanifolds where solutions behave differently. 
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Example 3.10. For F := C(x, y), A = {^, ^-} and U = {u} consider 

(xy - 1) • u(x,y) = 

and determine solutions on C 2 and its submanifolds. A differential THOMAS decomposition over 
F{u} simply reproduces this equation, because (xy — 1) G F \ {0}. However, we can model a 
search for solutions on submanifolds by adding two differential indeterminates X and Y to U 
and consider the equations. In order to allow splitting the manifold C 2 , we add two differential 
indeterminates X and Y to U which model the A- field elements x and y. Thus, we have to 
consider the additional equations (X x ■ (X x — 1)) = , (X y ) = , (Y y ■ (Y y — 1)) = , (Y x )= together with the 
modified equation ((XY — 1) • u) = . A differential THOMAS decomposition with X,Y <C u yields 
three systems: 

(i) (XY-1) = , (X x ) =1 (X y ) = , (X)? 

(ii) (u) =) (Y x ) = , (Yy(Y y -l)U (X X -(X X -1)U (X y ) = , (X)*, (XY-1)? 
(in) (u)=, (Y x ) = , (Y y ■ (Y y - 1)) = , (X) = 

System (i) allows an arbitrary function u on the submanifold M C C 2 defined by xy — 1 = as a 
solution. The other systems (ii) and (Hi) determine u = as the only solution on C 2 \ M . 

4. Implementation 

In this section, we describe our implementation of the decomposition algorithm. First, we 
list some other implementations of triangular decomposition algorithms. Second, we give some 
typical optimizations to make the computations feasible. Third, we describe our implementation 
in Maple. Fourth, we give benchmarks to get a more detailed and practical comparison between 
different decomposition algorithms. 

4-1. Implementations of Similar Decomposition Algorithms 

The RegularChains package by Lemaire et al. (2005) is shipped with recent versions of Maple. 
It contains the Triangularize command, which implements a decomposition of an algebraic variety 
given by a set of equations by means of regular chains. If the input also contains inequations, the 
resulting decomposition is represented by regular systems instead. It is possible to make these 
decompositions disjoint using the MakePairwiseDisjoint command. 

The epsilon package by Wang (2003) implements different kinds of triangular decompositions 
in Maple. It is the only software package besides our own that implements the algebraic Thomas 
decomposition. It closely resembles the approach that Thomas (1937, 1962) suggested, i.e., poly- 
nomials of higher leader are considered first. All polynomials of the same leader are combined into 
one common consequence, resulting in new conditions of lower leader. These are not taken into 
account right away and will be treated in later steps. Contrary to our approach, one cannot reduce 
modulo an unfinished system. Therefore, one needs extra inconsistency checks to avoid spending 
too much time on computations with inconsistent systems, epsilon implements such checks in order 
to achieve good performance. 

The Maple packages diffalg by Boulier and Hubert (1996-2004) and DifferentialAlgebra by 
Boulier and Cheb-Terrab deal with ordinary and partial differential equations as described by 
Boulier et al. (2009). They compute a radical decomposition of a differential ideal, i.e., a descrip- 
tion of the vanishing ideal of the Kolchin closure (Kolchin, 1973, §IV.l) of the set of solutions. 
Computation of integrability conditions is driven by reduction of A-polynomials (Rosenfeld, 1959, 
Sect. 2), which are the analogon of s-polynomials in differential algebra. Just like in RegularChains, 
this approach usually does not give disjoint solution sets, although in principle disjointness might 
be achieved. The diffalg package has been superseded by DifferentialAlgebra in Maple 14. Differ- 
entialAlgebra is based on the BLAD-libraries by Boulier (2004-2009) which have been designed as a 
set of stand-alone C-libraries with an emphasis on usability for non-mathematicians and extensive 
documentation. 
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4-.2. Algorithmic Optimizations 

In this subsection, we describe algorithmic optimizations helpful for a reasonably fast imple- 
mentation of the Decompose algorithm. 

In our algorithm, pseudo remainder sequences for the same pairs of polynomials are usually 
needed several times in different branches. As these calculations are expensive in general, our 
implementation always keeps the results in memory and reuses them when the same pseudo re- 
mainder sequence is requested again to avoid repeated computations. 

Coefficient growth is a common problem in elimination. Polynomials should be represented 
as compact as possible. Once we know that the initial of a polynomial is non-zero, the content 
of a polynomial (in the univariate sense) is non-zero, too. Thus, every time an initial is added 
to the system as an inequation, we can divide the polynomial by its content. Additionally, the 
multivariate content, which is an element of F, can be removed. 

The reduction algorithms (2.6) and (3.3) do not recognize that non-leading coefficients are 
zero. However, we can reduce the coefficients modulo the polynomials of lower leader, in addition 
to reduction of the polynomial itself. Thereby, in some cases the sizes of coefficients decrease, 
in other cases they increase. The latter is partly due to multiplying the whole polynomials with 
the initials of the reductors. Finding a good heuristic for this coefficient reduction is crucial for 
efficiency. 

Factorization of a polynomial improves computation time in many cases. More precisely, the 
system <SW{(p-g)=} decomposes disjointly into (SU{p = }, SU{p^, q=}) and the system S\J{(j)-q)jt} 
is equivalent to S U {p^,q^}. In most cases, the computation of two smaller problems resulting 
from a factorization is cheaper than the computation of the big, original problem. This idea 
extends to factorizations over an extension of the base field: Let Yi := |xj | Xj < a;,, {St)x } 7^ $| 

and Zi := |xj | Xj < Xi, (St)^ = Assume that (Sr)^ is irreducible over the field Fj := 

F(Zi)[Yi]/ ((ST)< Xi ) for all i £ {1, ...,n.}, where ((Sx)^) i s the ideal generated by (St)< X( in 
the polynomial ring F(Zi)[Yi]. Factorization over F n instead of F may split the polynomial into 
more factors, but it is not clear whether this improves runtime. Preliminary tests show that 
factorization over F should be preferred for F = Q. 

In the algebraic algorithm, polynomials need not be square- free when they are inserted into the 
candidate simple system. Efficiency can sometimes be improved by postponing the computation 
of the square-free split as long as possible. However, this is not possible for the differential case. 
Differential polynomials need to be made square-free to ensure that their separant is non-zero, i.e. 
non-trivial prolongations have a non-zero initial. 

In the differential case, application of criteria can decrease computation time by avoiding 
useless reductions of non-reductive prolongations. Janet's combinatorial approach already avoids 
many reductions of A-polynomials, as used in other approaches (see Gerdt and Yanovich (2006)). 
In addition, we use the involutive criteria 2-4 (cf. Gerdt and Blinkov (1998a); Gerdt (2005); 
Apel and Hemmecke (2005)), which together are equivalent to the chain criterion. Applicability 
of this criterion in the non-linear differential case was shown in (Boulier et al., 2009, §4, Prop. 5). 

The axioms of a selection strategy (see Definition (2.22)) already strongly limit the choice 
for the polynomial considered in the current step. However, the remaining freedom is another 
important aspect for the speed of an actual implementation. We will describe different selection 
strategies in §4.3 and compare them in the benchmarks. 

As described up to now, the algorithm often keeps on computing with inconsistent systems. We 
want to optimize the algorithm to detect the inconsistencies as early as possible. This allows the 
algorithm to discard inconsistent systems as early as possible. One of the problems are selection 
strategies that postpone the costly treatment of inequations. A test to detect whether inequations 
in Sq reduce to zero is comparably cheap. 

Another possible improvement is parallelization, since the main loop in Decompose (2.25) can 
naturally be used in parallel for different systems. 

4-. 3. Selection strategies 

We consider our two main approaches to selection strategies (see Definition (2.22)). 
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1. The "equations first" strategies: Select only chooses an inequation if Q does not contain any 
equations. Among the equations or inequations, it prefers the ones with smallest leader. 

2. The "leader first" strategies: Select always chooses an equation or inequation with the small- 
est leader occurring in Q. If there are both equations and inequations with that leader, it 
chooses an equation. 

In both approaches, if the above criteria do not yield a unique choice, we compare the leader of 
the initial and choose the smaller one. We apply the last test recursively to the initial of the initial 
and so on. At this point, it is still possible that we fail to make a unique selection. However, these 
cases are rare and there does not seem to be a considerable performance advantage for any choice. 
Therefore, it suffices to make an arbitrary (but preferably unique) choice. 

In our experimental observation "leader first" strategies usually produce decompositions with 
less systems, while "equations first" strategies are more efficient (cf. §4.5). 

4-. 4- Implementation in MAPLE 

Both the algebraic and the differential case of the Thomas decomposition algorithm have been 
implemented in the Maple computer algebra system. Packages can be downloaded from our web 
page (Bachler and Lange-Hegermann (2008-2011)), documentation and example worksheets are 
available there. 

The main reason for choosing Maple for the implementation is the collection of solvers for 
polynomial equations, ODEs, and PDEs already present. Furthermore, fast algorithms exist for 
polynomial factorization over finitely generated field extensions of Q and for gcd computation. 

The AlgebraicThomas package includes procedures to compute a Thomas decomposition, re- 
duce polynomials modulo simple systems and compute counting polynomials (cf. Plesken (2009a)). 
Furthermore, it can represent the complement and intersection of solution sets as decompositions 
into simple systems. Finally, a comprehensive Thomas decomposition can be computed, this topic 
will be discussed in a later publication. 

Example 4.1. We demonstrate how to use the AlgebraicThomas package by computing a decom- 
position of the system in example (2.5). 

> wzth( AlgebraicThomas) : 

> p := a*x~2 + b*x + c; 

p := x 2 a + x b + c 

> S := AlgebraicThomasDecomposition( [p] , [x, c,b , a] ) ; 

S := [[x 2 a + xb + c = 0, 4 ca - b 2 ^ 0, a ^ 0], [2x a + b = 0, 4 ca - b 2 = 0, a ^ 0], 

[xb + c = 0, b^0, a = 0], [c = 0, b = 0, a = 0]] 
Information about leader and main degree can optionally be included in the output. 

> map (printSy stem, S, ["PT", "LR "] ) ; 

[[[x 2 a + xb + c = 0, x 2 }, [ica-b 2 ^ 0, c], [a ^ 0, a]], 
[[2xa + b = 0, x], [Aca-b 2 =0, c], [a ^ 0, a]], 

[[xb + c = 0, x], [6^ 0, b], [a = 0, a]], [[c = 0, c], [6 = 0, b], [a = 0, a]]] 
It is possible to include inequations in the input to exclude some degenerate cases: 

> q := a<>0; 

q:=a^0 

> T := AlgebraicThomasDecomposition( [p, q] , [x,c,b ,a] ) ; 

T := [[x 2 a + xb + c = 0, 4 ca - b 2 ^ 0, a ^ 0], [2 x a + b = 0, 4 ca - b 2 = 0, a ^ 0]] 

Features for the differential package DifferentialThomas include arbitrary differential rankings, 
using special functions implemented in Maple as differential field elements, computation of power 
series solutions, and a direct connection to the solvers of Maple for differential equations. 
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Example 4.2. We treat the following control theoretic example taken from Diop (1992). 

> with(DifferentialThomas): 

> ComputeRanking ( [t] , [x2,xl ,y,u] , "EliminateFunction") ; 

This creates the differential polynomial ring Q{x^ 2 \ x^\ y, u} for A = {■§[}■ Here u indicates 
the input, x^- 1 ' and x^ the state, and y the output of the system. The chosen ranking "<" is the 
elimination ranking with x^ 3> a;' 1 ^ > y > i.e., x^ > x^ > j/k > u\ for all i,j,k, 1 G Z>o- 

> L:=[xl[l]-u[0]*x2[0] ,x2 [1] -xl [0] -u[0] *x2 [0] ,y[0]-xl[0]] : 

We follow (Diop, 1992, Ex. 1) and compute the external trajectories of a differential ideal 
generated by L, i.e. intersect this differential ideal with Q{y,u}. 

> res:=DifferentialThomasDecomposition(L, []) ; 

res := [DifferentialSystem, DifferentialSystem] 

We show the equations and inequations of the differential systems not involving x^ 1 ' or x^ . 
The chosen ranking guarantees that the systems shown determine the external trajectories of the 
system: 

> PrettyPrintDifferentialSy stem (res [1] ) : remove (a->has (a, [xl , x2] ) ,%) ; 

\-<t) y(t)) + a y(*)) + a y (t)) a u®) + y w u^ 2 = o, + o] 

> PrettyPrintDifferentialSy stem (res [2] ) : remove (a->has (a, [xl , x2] ),'/,) ; 

[&y(t) = Q, u(t) = 0] 

These systems, having disjoint solution sets, are identical to the ones found in Diop (1992). 
4-5. Benchmarks 

In this subsection, we compare our two Maple packages to the other implementations of 
triangular decompositions mentioned in §4.1 using benchmarks. Not all of the implementations 
compute equivalent results. This should be considered when comparing the timings. We omitted 
examples where all tested systems took less than one second to complete the computation or could 
not be computed by any software package. 

All benchmarks have been performed with Linux x86-64 running on a third generation Opteron, 
2.3 GHz. The time limit has been set to 3 hours and available memory is limited to 4 GB. All times 
are given in seconds. The polynomial multiplication in Maple 14 benefits from a new parallel 
implementation (cf. Monagan and Pearce (2009)). Nonetheless, we state the total CPU time in 
our benchmarks, as returned by Maple's time command. 

By default, both of our Maple packages behave as follows: 

• Polynomials are factorized over Q. 

• The content of polynomials is removed. 

• The selection strategy is an "equations first" strategy, as decribed in §4.3. 

• After reducing a polynomial, we always reduce its coefficients fully. 

• Inequations in Sq are reduced for early inconsistency checks. 
See §4.2 for details. 
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Table 1: Comparison of algebraic decompositions 1: polsys50 from Wang (2003) 
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4-5.1. Algebraic Systems 

For testing the AlgebraicThomas package, we used two sets of examples, namely, the test ex- 
amples from the polsys50 file in Wang's epsilon package (Wang (2003)) printed in table 1 and 
the examples from Chen et al. (2007) as shown in table 2. 

In contrast to Algorithm 2.25, the implementation in the AlgebraicThomas package inserts 
equations or inequations into St without making them square- free first. It delays this computation 
as long as possible, sometimes until the end of the decomposition. This avoids some expensive 
and unnecessary discriminant computations entirely. 

We compared AlgebraicThomas with the RegularChains package from Maple 14 and epsilon. 
We also tested the AlgebraicThomas and RegularChains packages in different configurations. The 
timings in Maple 14 of the following procedures are being compared: 

• (RC1) RegularChains [Triangularize] . 

• (RC2) RegularChains [Triangularize] with the 'output ' = 'lazard' option set. 

• (RC3) RegularChains [Triangularize] with the 'radical ' = 'yes' option set. 

• (DW1) epsilon [RegSer] . 

• (DW2) sisys [simser] . 

• (ATI) AlgebraicThomasDecomposition. 

• (AT2) AlgebraicThomasDecomposition with factorization disabled. 

• (AT3) AlgebraicThomasDecomposition with a "leader first" selection strategy (cf. §4.3). 
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Table 2: Comparison of algebraic decompositions 2: Test examples from Chen et al. (2007) 



Name 








JJW1 


JJWz 


ATI 
All 


Al A 


Al 


Al 4 


AlkashiSinus 


0.6 


0.6 


0.8 


0.1 


7.1 


5.7 


2.6 


6.5 


3.6 


Bronstein 


0.4 


0.5 


0.5 


0.2 


0.4 


0.3 


0.4 


1.1 


0.4 


Cheaters-homotopy-easy 


0.7 


> Sh 


532.5 


> 4GB 


> 4GB 


> Sh 


> 4GB 


> 3h 


> 3/i 


Cheaters-homotopy-hard 


0.7 


> Sh 


559.8 


> 4GB 


> 4GB 


> Sh 


> 4GB 


> Sh 


> Sh 


Gerdt 


1.4 


1.4 


1.4 


2.0 


2.2 


8.1 


3.2 


0.5 


1532.1 


Hereman-2 


0.8 


1.0 


0.8 


0.3 


0.4 


0.3 


1.2 


0.3 


0.5 


Hereman-8-8 


26.9 


31.6 


208.3 


> Sh 


> 3h 


> 3/i 


> Sh 


> 3/i 


> 4GB 


rvd V 


700 O 

1 zz.z 


1C\1 1 

f U / .1 


/ ZD. 1 


> on 


> OIL 


> Oil 


> Oil 


> Oil 


> Oil 


Lanconelli 


0.4 


0.6 


0.4 


0.2 


0.4 


0.4 


1.3 


0.3 


0.3 


Lazard-ascm2001 


1 2 


1 l .0 


1 4 


s> Oil 




\ AC ft 




s> Oil 


^> OIL 


Ley kin- 1 


5.6 


8.0 


5.8 


> 3h 


> 3h 


2.3 


> 3h 


6.0 


1.4 


Maclane 


2.4 


6.7 


2.6 


3576.5 


> 4GB 


7.4 


17.4 


13.1 


7.0 


MontesSlO 


0.5 


1.0 


0.6 


0.4 


17.7 


2.0 


2.2 


2.3 


1.5 


MontesSll 


0.2 


0.5 


0.2 


0.2 


> 3h 


23.5 


> 3h 


21.9 


12.4 


MontesS12 


0.5 


3.2 


0.5 


1.6 


> 3h 


9.4 


31.0 


13.6 


115.4 


MontesS13 


0.3 


0.5 


0.3 


0.2 


0.5 


0.8 


1.8 


1.1 


0.9 


MontesS14 


0.7 


1.3 


0.8 


> 3h 


> 3h 


6.0 


> 4GB 


14.5 


12.1 


MontesS15 


1.2 


1.8 


1.3 


0.6 


8.2 


4.8 


3.8 


6.7 


3.3 


MontesS16 


4.3 


3.4 


4.2 


1.5 


1.5 


2.5 


2.2 


3.2 


1.5 


MontesS7 


0.4 


0.5 


0.4 


0.2 


0.5 


0.7 


0.6 


2.5 


1.2 


Neural 


0.5 


0.7 


0.6 


> 3h 


> 4GB 


1.4 


3050.9 


1.7 


1.2 


Pavelle 


1.1 


15.9 


1.4 


> 3h 


> 4GB 


> 3h 


> 3h 


> 3h 


> 3h 


Wang93 


1.2 


1.3 


1.2 


1.6 


3.4 


4.9 


> 3h 


3.4 


6.5 


genLinSyst-3-2 


0.3 


1.1 


0.3 


0.2 


0.2 


0.3 


0.2 


0.3 


0.2 


genLinSyst-3-3 


0.3 


4.5 


0.4 


1.2 


1.2 


6.0 


2.5 


4.8 


1.1 



• (AT4) AlgebraicThomasDecomposition with coefficient reduction disabied. 

We compare table 1 and 2 within our own implementation. We observe, that (AT2) is much 
slower than (ATI) and, thus, conclude that factorization is vital to make many computations 
feasible. In a few examples, we see the relative advantage of the default selection strategy compared 
to the one used in (AT3). Generally speaking, disabling coefficient reduction increases computation 
time for (AT4) , but there are some strong counterexamples to this observation. This indicates that 
different strategies for coefficient reduction, as seen in (ATI) and (AT4), should be investigated 
further. 

The programs sisys [simser] (DW2) and AlgebraicThomasDecomposition (AT1-4) are the 
only ones that compute a Thomas decomposition. All test examples that could be computed by 
(DW2) could also be computed by (ATI). However, there are some examples that RegularChains 
(RC1) or epsilon (DW1) could treat, but we could not decompose into simple systems. Moreover, 
the test examples indicate that (RC1) is in general faster than (ATI) in the positive-dimensional 
case. Our evaluation suggests that this is due to the strict square-free property of simple systems. 
In the zero-dimensional case, however, the situation is less clear, since there are examples where 
(RC1) is faster than (ATI) and vice versa. 

4-5.2. Differential Systems 

We compared DifferentialThomas with the packages diffalg and DifferentialAlgebra. Finding a 
suitable set of benchmark examples for the differential case was more difficult. We are not aware of 
any sets of standard benchmarks. Thus, we used a collection of examples, which we came across 
in our work. These examples are published on our homepage (Bachler and Lange-Hegermann 
(2008-2011)). 

The timings of the following procedures are being compared: 

• (DA) Dif f erentialAlgebra [Rosenf eld_Groebner] . 
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Table 3: Benchmarks for ODE systems 



Name 


DA 


da 


DTI 


DT2 


DT3 


DT4 


Diffalg4 


2.9 


2.9 


852.5 


> 3h 


8932.4 


36.0 


LLG3 


0.5 


> 3h 


5.4 


5.6 


4.4 


4.9 


LLG4 


0.3 


19.1 


2.6 


37.4 


20.3 


4.0 


ODE1 


2.4 


3.7 


0.6 


0.3 


0.6 


0.8 


ODE6 


2.3 


1.5 


0.8 


1.2 


0.6 


0.8 


ODE7 


> 3h 


> 3h 


3.2 


60.8 


47.2 


5.4 


kepler vs newton 


0.7 


0.9 


1.4 


2.8 


0.8 


2.0 


kepplerl 


0.1 


0.2 


1.1 


0.6 


0.8 


1.1 


keppler2 


0.1 


0.1 


1.3 


1.5 


1.1 


0.8 


keppler3 


0.1 


0.2 


1.1 


0.6 


0.7 


1.1 


murrayl 


0.1 


0.4 


1.9 


2.4 


1.7 


2.2 


murray2 


0.1 


0.1 


0.7 


1.4 


0.8 


0.6 



Table 4: Benchmarks for PDE systems 



Name 


DA 


da 


DTI 


DT2 


DT3 


DT4 


Cyclic 5 variantl 


> 3h 


2.9 


1.3 


1.5 


1.4 


1.2 


Cyclic 5 variant2 


> 3h 


> 3h 


2.3 


2.6 


0.7 


2.3 


Diffalg2 


0.5 


0.3 


1.4 


1.5 


41.5 


1.2 


Diffalg3 


0.2 


0.5 


0.9 


0.6 


1.6 


0.7 


Ibragimov 2, 17.9c 


> 3h 


> 3h 


18.4 


> 3h 


40.8 


12.3 


PDE6 


> 4GB 


> 4GB 


11.1 


23.0 


> 4GB 


16.5 


PDE7 


> 4GB 


116.7 


91.3 


83.3 


> 4GB 


41.4 


PDE8 


> 4GB 


> 4GB 


6.8 


> 4GB 


14.2 


7.5 


Riquier lb 


0.1 


0.1 


1.9 


1.9 


1.9 


2.0 


Riquier 3a 


0.3 


0.2 


0.8 


1.1 


0.5 


0.6 


Riquier 3b 


0.6 


0.6 


1.4 


1.7 


1.4 


1.2 


boulier 


> 3h 


546.4 


2.5 


1690.9 


2.6 


1.5 


cyclic 6 


> 4GB 


571.9 


160.7 


159.8 


349.6 


154.1 


noon6 


> 4GB 


72.6 


40.6 


36.8 


63.7 


31.0 



• (da) dif f alg [Rosenf eld_Groebner] 12 . 

• (DTI) Dif f erentialThomasDecomposition. 

• (DT2) Dif f erentialThomasDecomposition with factorization disabled. 

• (DT3) Dif f erentialThomasDecomposition with a "leader first" selection strategy (cf. §4.3) 

• (DT4) Dif f erentialThomasDecomposition with coefficient reduction disabled. 

We want to mention one further example not included in the benchmark table, ft is the test 
example 5 of dif f alg. None of the packages in their default setting could compute this example. 
Still, diffalg and Dif f erentialAlgebra were able to do so instantaneously by a change of 
ordering strategy. 

The comparison between (DTI), (DT2) and (DT3) is similar to the algebraic case. In partic- 
ular, factorization should be enabled and the default selection strategy should be preferred. In 
contrast to the algebraic implementation, the comparison of (DTI) and (DT4) is less conclusive. 

All test examples which could be computed by Dif f erentialAlgebra or diffalg could also 
be computed by our default strategy (DTI). For ODEs, the three packages show similar timings, 
but for PDEs, Dif f erentialThomasDecomposition appears to be faster. This might be explained 
by the involutive approach, which we utilize to make the subsystems coherent. A similar result 
can be found for the GINV-project (cf. Blinkov et al. (2010a), Blinkov et al. (2010b)). 



with _env_dif f alg_uses_Dif f erentialAlgebra: =f alse 
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